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1. Introduction. As a general rule, a problem in physics is not considered solved 
unless the solution can be expressed in analytical form. The same usually holds true 
in the case of engineering problems, although there the art often progresses faster than 
the theory under the impact of economic forces, and the engineer is often forced to seek 
a solution by means of an experimental setup, or possibly by means of some numerical 
or graphical process. ; 

The disadvantage of a numerical or graphical method is its lack of generality, its 
tendency towards inaccuracy, particularly owing to cumulative errors, and its in- 
ability to exhibit optimum values for the parameters involved, particularly if these 
have to be in numerical rather than in symbolic form. On the other hand, these meth- 
ods often yield answers to problems that the analytical method cannot handle, and 
furthermore are often very effective as teaching aids. This is particularly true of the 
graphical methods. 

It is the purpose of this article to illustrate the application of graphical construc- 
tions to problems involving nonlinear circuits, particularly those containing vacuum 
tubes. It is the writer’s hope that some mathematician will be sufficiently attracted 
to this method to attempt to establish it on a more general basis, possibly something 
akin to the collection of theorems of ordinary Euclidean or of Projective Geometry. 

2. Definition of graphical method. Before proceeding with a description of the 
method it will be desirable to define it. By graphical constructions are meant those 
geometrical manipulations by which a solution to a problem is obtained. It may be 
necessary to slide a curve representing a relationship between two variables along 
the axis of the independent variable, and to find (geometrically) where it intersects 
another curve representing a second relationship between the two variables. The 
manipulations may be more involved than those of simple translation along the axis, 
and it is to be stressed that the restriction of ruler and compass constructions is not 
invoked in these manipulations. 

It is apparent that the method is not that usually understood by the average engi- 
neer, namely, the plotting of a complicated generalized analytical expression to per- 
mit values to be taken off the graph in order to obviate the need for computing the 
value of the expression every time the problem arises. 

3. Simple series nonlinear circuit. As an elementary example of a graphical con- 
struction, let us consider the circuit shown in Fig. 1, that of a diode (two-element 





* Received March 13, 1945. 








186 ALBERT PREISMAN [Vol. III, No. 3 


vacuum tube) in series with a resistance R and a source of d.c. potential E. It is de- 


sired to find the current flow in this circuit. 

It is necessary to know the voltage-current e—7z relationship for the diode, and 
for the resistor R. We assume for simplicity that the latter is a linear resistance. 
Then the e—i relationship is that shown in Fig. 2. The curve is a straight line making 
an angle @ with the voltage axis, such that 

cot 6 = R, (1) 


the resistance of the device. This slope is constant and hence R has a fixed value: 


so many volts per ampere, or ohms. 
On the other hand, the diode has the characteristic shown in Fig. 3. Here, for 
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negative values of voltage (piate negative to cathode) no current can flow; while for 
positive values of voltage current flows in such manner as to generate the curve shown. 
The ideal diode would have the following equation for positive plate voltages 

i = Aed/?, (2) 
but actual diodes depart to some extent from the above equation owing to such fac- 
tors as initial velocity of emission of electrons from the cathode, the effect of the sup- 
porting members for the cathode and plate, etc. 

The diode is a nonlinear device; first because of the break in the curve at the origin 
and second because even for positive plate voltages the e—7 relationship is usually 
not a straight line. One can define the resistance as 

1) the reciprocal slope of the secant line to any point of the curve (this is the 
so-called d.c. resistance) or 

2) the reciprocal slope of the tangent line to any point of the curve (this is 
usually called the a.c., incremental, or variational resistance of the device). 
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Such concepts have limited utility however, since the resistance in either case is no 
longer a constant, but a function of the applied voltage or current through the device. 
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The graphical method to be described takes the fundamental e—i relationship, or 
terminal characteristic as it has been called by Kirschstein, and operates directly with 
it. Furthermore, the curve does not have to be analytic, nor even expressible in the 
form of an equation; it can simply be a plot of experimental data, although this in- 
volves interpolation between experimentally determined points. 

The process of finding the current through the two in series for the impressed 
voltage E is essentially that of solving the two equations for the terminal character- 
istics simultaneously under the condition that the sum of the voltage drops across 
the two elements must equal the impressed voltage EZ. Thus, if the relationship for 
the one element is i=fi(e), then that for the other is 1=f.(E—e), and it is desired to 
find a common value of i that satisfies both relationships. 

Since one or both of the above equations may be of degree higher than unity, the 
analytical solution cannot be effected by the method of determinants, but rather by 
the method of substitution, and finally results in the necessity for solving an equation 
of degree higher than unity. 

This, however, assumes that terminal characteristics can be represented by power 
series. The graphical method requires no such condition; it operates on the graphical 
plots directly. Thus, suppose the terminal characteristic of the diode is represented 
by AOB, Fig. 4. Let OC represent the magnitude of the impressed voltage E. Through 
C draw DC at an angle 0, as shown, such that cot @=R. Then the intersection of CD 
and AOB in D represents the required solution, in that DF is the common current 
in this series circuit; OF is the voltage drop in the diode; FC is the voltage drop in 
the resistor R; and clearly OF +FC equals the impressed voltage E. If E varies with 
time, DC can be shifted back and forth along the voltage axis at positions correspond- 
ing to the instantaneous values of E, and the intersections will furnish the corre- 
sponding instantaneous values of the current. 

The above solution represents a well-known method for solving two equations 
simultaneously when the equations are of degree higher than the first or even of 
transcendental nature. It will be of interest, however, to see how this method is ap- 
plied to a more complicated circuit. 

4. Triode tube and resistance in series. The next example will be that of a three- 
element or triode tube in se- 
ries with a resistance and a Triode 
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The resulting plot is a curved surface in space. It can be represented in two di- 
mensions by a family of curves which represent discrete projections of this surface 
upon any one of the three coordinate planes. For the problem at hand the most useful 
set of projections is that upon the e,—i, coordinate plane, in the form of a family of 
€p—ip curves with e, as the parameter. This is shown in Fig. 6 (solid lines). Curves 
for which e, is positive have been omitted for simplicity. 

Assume further for simplicity that Rz, the load resistance, is linear. The current 
through it is a function of but one voltage, that which must be applied across its 
terminals to produce the above current flow. To represent its terminal characteristic 

in three dimensions, it is plotted as a 


: & plane whose intersection with the e,—é, 
Up hy ¢ coordinate plane is a straight line paral- 
o lel to the e, axis. In this way the current 

A hy : 


in it is independent of the e, coordinate, 








e and is a (linear) function of but one 
voltage, that corresponding to the plate 
0 voltage e, of the tube. All points of this 
O c \ £s, & plane representing Ry, project over to 
a the 1, —e, coordinate plane as a straight 
“IG. 0. 


line that is also the intersection of the 
above R, plane with the i,—e, plane. 
The straight line makes an angle @ with the e, axis such that cot 0=R,, i.e., the 
R, plane is inclined at the angle @ to the e,—e, coordinate plane. 

The graphical solution consists in drawing the line of intersection EA at the angle 
6 to the e, axis. The intersection of EA with the tube family of curves gives the com- 
mon value of current flowing through the plate circuit of the triode and R, in series, 
for any given value of grid voltage e,. For example, at a moment when the signal 
voltage e, is passing through zero, the instantaneous value of the grid voltage e, is 
simply that of the bias battery, E.. The instantaneous value of the plate current is 
BC, where B is the intersection of AC with that curve of the plate family for which 
e,=£,. It is further to be noted that the instantaneous plate voltage e, is OC, and 
the instantaneous value of the voltage drop across R, is EC. 

For other instantaneous values of e,, other curves of the plate family are involved, 
and the process of determining the instantaneous values of plate current, plate volt- 
age, and load voltage (across R,) is identical to that described above. Thus, for a 
signal voltage impressed upon the input or grid circuit, the output signal voltage be- 
tween the plate and ground can be found. Such matters as the amplification of the 
stage, distortion in the output, etc., can then be determined. 

In passing, we may note here that the locus of the plate current for various values 
of e, is the intersection of the tube surface and the R, plane in space. This intersection 
is a curve in space, but fortunately its projection on the e,—7, plane is a straight line, 
namely the intersection of the R, plane itself with the e, —7, plane. It is for that rea- 
son that the e,—7, family of the tube curves is employed; the graphical construction 
is simply the points of intersection of a straight line representing Rz, with the above 
plate family. 

The above problem can become much more complicated under certain conditions. 
For example, if the input signal voltage is great enough, the grid can be driven posi- 
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tive with respect to the cathode, whereupon it draws current during the positive peak 
of the a.c. cycle. If the signal source has appreciable internal impedance, then a volt- 
age drop will occur in the source during the above portion of the cycle, and the actual 
voltage applied to the grid will differ from the generated voltage e,. 

It is therefore necessary to determine the actual grid voltage before the plate 
current can be found. Another complication arises however, in that the grid current 
(and hence the actual grid voltage) is a function not only of the positive grid voltage, 
but of the plate voltage as well. This is because the space current divides between 
the two electrodes in a manner depending upon the two electrode voltages. At the 
same time the plate voltage is a function of R, and the grid voltage. Thus the above 
simple graphical construction can become quite involved if merely the input signal is 
increased to a point where the grid is driven positive. 

5. The balanced amplifier. Instead of investigating such details, important though 
they may be, it will be of interest to examine another type of circuit very important 
in the communication industry. Reference is made to the push-pull or balanced am- 
plifier. The circuit is shown in Fig. 7. 

In (A) is shown the actual circuit, whereas in (B) is shown an idealization or 
equivalent form better suited for the purpose of analysis. In the actual circuit (A), 
two tubes are employed, inductively 
coupled to each other and the output /put, 
load resistance r;, by an output trans- 
former. The signal on one grid is 180 de- 
grees out of phase with that on the other 
grid, as is suggested by the symbols +e, 
and —e,. The bias voltage E,, on the other 
hand, is applied to both grids in the same 
polarity; and the plate supply voltage is 
applied to the two tubes in the same polar- A 
ity too, as shown. 

The actual load resistance rz, and the ——> J, 
output transformer can be replaced by the Lith said 
center-tapped inductance and reflected 
load resistance R, as far as the tubes are a | 
concerned. The simplified circuit is shown jaleloje os zy 
in (B), Fig. 7. In using this equivalent cir- 
cuit, it is tacitly assumed that the actual a 
output transformer is an ideal transformer —> I, 
having infinite primary and secondary 
open-circuit inductance, no distributed ca- 
pacity, unity coefficient of coupling be- Fic. 7. 
tween windings, etc. In the equivalent 
circuit the center-tapped inductance is assumed to be infinite in value and to have 
unity coupling between the two halves of the complete winding. Ordinarily this is a 
reasonable assumption. 

As a result, the current in one-half of the winding cannot at any moment exceed 
that in the other half for otherwise an infinite counter-electromotive force would be 
induced in the windings that would tend to prevent such an unequality from taking 
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place. The currents in the windings can vary, however, provided they remain equal 
to one another at all times. Finally, two further assumptions are made, namely, that 
the signal voltages e, and —e, applied to the two grids are at all times equal and op- 
posite to one another, and that the two tubes have identical terminal characteristics. 
These two assumptions seem also reasonable. 

Consider first that e, equals zero (no signal is applied). The bias on each grid is E,, 
and the plate voltage for either tube is Ey», hence the two plate currents J; and J; are 
equal to one another. Since they flow in opposite directions from the ends of the wind- 
ing to the center tap, they balance each other magnetically in the output inductance 
and produce no voltage across the ends. Consequently no current flows in the load 
resistance R;. 

Now suppose that a signal voltage is impressed such that the top grid is driven 
positive by an amount e,; from its normal d.c. negative bias value of E., and that the 
bottom grid is driven more negative by an equal amount, i.e., —é,1. The two plate 
currents will now vary in opposite directions, namely, J; will increase and J, will de- 
crease. However, the sum of these two currents flows through the plate power supply, 
and owing to the infinite inductance of the center-tapped winding, (J:+J2)/2 flows 
down through the top half of the winding, and an equal amount flows up through the 
bottom half, to combine at the center tap to furnish the sum (1,;+J2) flowing through 
the power supply. 

Since (J:+J2)/2 is the average between J; and J», it is equal to neither, and from 
the principle of continuity of current flow, the difference 


hh- 31+ 12) = (1 + 12) — Tn = 3(1 — 12) (3) 


must flow through R;. A quick check will indicate that Kirchhoff’s current law is 
satisfied at each junction. 

The current (J;—J2)/2 is the output current. In flowing through R_,, it sets up a 
voltage drop 

Ex = 3(11 — I2)Rt. (4) 

Half of this or E,/2 appears across each half of the output winding of such polarity 
that the instantaneous plate voltage of the top tube is E,—(£,/2) and that of the 
bottom tube is E,+(E,/2). 

Thus the following facts have been brought to light: 

1) The grid voltages change by equal but opposite increments from their 
common bias value E, owing to the center tap on the input transformer secondary. 

2) The plate voltages change by equal but opposite increments from their 
common supply value E,, owing to the center tap on the output inductance. More- 
over, the plate voltage increments are opposite in sign to the corresponding grid 
voltage increments. 

3) The plate currents change in opposite directions in the same sense as the 
corresponding grid voltages, but not necessarily to an equal degree. If the tubes 
are nonlinear, as is usually the case, then the increase in plate current of either 
tube for a positive increment in grid voltage is not necessarily the same as the de- 
crease in plate current for an equal negative increment in grid voltage. 

From the above facts several graphical constructions are available to determine 
the plate current and plate voltage variations in the tubes, the output current and 
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voltage, the power output, and the d.c. power input. The following graphical method 
is preferred by the author. In Fig. 8 is shown the plate family of curves for either 
tube. If there is no signal input, the only voltages present are the d.c. potentials En 
applied to the two plates and E, applied to the two grids. The current through either 
tube is then 7, =E p»B, a direct current. 


Up 











Now suppose that equal and opposite signal voltages e, and —e, are applied to the 
grids in addition to E,. Then the current in the one tube will increase from BEp» 
to DG, and that in the other tube will drop to FH, as shown. The plate voltage of the 
first tube will drop from OE,» to OG = (Ey —Ae,), and that in the other tube will rise 
by an equal amount to OH = (Ew +Ae,). 

It is also clear from Fig. 8 that DJ represents the difference between the two cur- 
rents or (J;—J2), and JF represents 2Ae,, the voltage across the output inductance 
and previously denoted by E, in Fig. 7. From Eq. (4), it is evident that 


JF/DJ = E1/(h — 12) = R1/2. (5) 
Thus DF makes the angle @ with the e, axis such that 
cot @ = R,/2. (6) 


It is also evident from the geometry of the figure that DC =CF, i.e., that the ordinate 
through E}, bisects line DF in C. 

The above facts suggest the following method of graphical construction. We hold 
a rule at the angle @ and slide it up or down until the segment between the desired 
€, —i» curves (corresponding to equal and opposite grid voltage excursions from the 
bias value E,) is bisected by the ordinate through Ey». The intersections of the rule 
with the two e, —i, curves gives the two instantaneous values of the two tube currents 
I, and I», corresponding to the signal voltages e, and —e, and to the plate load re- 
sistance R,, or rather to R,/2. 

Then another pair of equal and opposite grid signal voltages are chosen, and the 
process repeated. This is continued until as many pairs of instantaneous grid signal 
voltages have been used as is desired. For a symmetrical signal voltage, such as a 
sine wave, instantaneous values for only one-quarter of a cycle are required. 

When the above graphical construction is performed, there is obtained a curve 
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on the plate family of curves such as that shown in broken lines ABCDE in Fig. 9. 
This represents the locus of the current for either tube over a cycle of grid signal volt- 
age. It also represents the terminal characteristic for R, as it appears to either tube 
in the presence of the other tube. 

The significance of the last statement is as follows: the two tubes may be regarded 
as two generators connected to a common load R,. Owing to their nonlinear charac- 
teristics, the tubes do not share the load equally throughout the signal cycle; that 
tube whose apparent internal resistance is lower takes a greater share of the load, 
i.e., furnishes more than half of the load current (J;—J2)/2 flowing through R,. As a 
result, R, appears as a variable or nonlinear resistance to either tube even though it 
is actually a linear resistance, and its terminal characteristic on either tube’s e,—7, 
family of curves is in itself a curved rather than a straight line. 


lp y 








~) 
» 





FIG. 9. 


Lack of space precludes a detailed discussion of this interesting circuit. However, 
several important features will be presented. As indicated in Fig. 9, the two e,—1p 
curves passing through B and D, respectively, represent equal and opposite grid 
swings. The corresponding currents J; and J, for the two tubes are BF and zero; in 
short, the tube experiencing the negative grid swing has just reached plate current 
cutoff. 

For ¢,—%, curves passing through A and E, corresponding to a still greater grid 
swing for either tube, J; is AG, and J, still remains zero. This means that the second 
tube is inoperative over this part of the cycle and acts therefore as if it were discon- 
nected. Under these conditions R; appears to the operative tube as R,/4, which can 
be expected since the 2 to 1 turns ratio of the output inductance will produce this 4 
to 1 impedance transformation if it is unhampered by the other tube. 

Portion BA is therefore a straight line whose reciprocal slope corresponds to R,/4. 
It is easy to show that if it were prolonged, it would pass through Ey». Normally 
the tubes are operated so that maximum grid signal voltage drives each tube alter- 
nately to cutoff or beyond. Maximum output occurs if R,/4 equals either tube’s 
apparent internal plate resistance at the peak of the cycle. The plate resistance of 
either tube is given by the reciprocal slope of the e,—1, curve at point A. Hence a 
quick determination for the optimum value of R,, or rather R,/4, is to draw a line 
through E;, at an angle equal to that of the e,—i, curve at point A, and calculate 
from the reciprocal slope of this line the value of R,/4 and hence of R,. The complete 
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characteristic can then be determined by means of the sliding rule as described pre- 
viously. 

Fig. 8 also reveals an interesting point. CE), is the average between DG and FH, 
i.e., it represents (J,;+J,2)/2. This is the mid-branch current that flows through the 
plate supply, as indicated in Fig. 7(B). For various pairs of values of J; and J: as 
determined by the sliding rule, the average, or (J1:+J2)/2 moves up and down along 
EppA. This is a vertical line or ordinate, and indicates that the resistance to the mid- 
branch current is zero. This has been tacitly assumed; the output inductance and the 
plate supply have been assumed to be free of resistance. 

If this is not the case, then a line must be drawn through E,» whose reciprocal 
slope indicates one-half the value of the mid-branch resistance that is present, and 
the sliding rule must be bisected by this line rather than the ordinate Ey»A, as is the 
case in Fig. 8. From this follows several further interesting characteristics.’ 

Another point is that not only is the locus of the mid-branch current along the 
ordinate E,»A in Fig. 8, but that this current executes two alternations per cycle of 
the grid signal voltage. This means that the mid-branch current is at least double the 
frequency of the incoming signal; actually, for perfect symmetry, all the even har- 
monics generated by the tubes flow in parallel through the mid-branch portions of 
the circuit, while the odd harmonics, including of course the fundamental, flow 
through the output resistance R,. Thus, if the tube characteristics are such that the 
second harmonic is quite prominent, but the third (and higher) harmonics are of small 
amplitude, then the output wave will be a fairly faithful copy of the input grid signal 
voltage and the stage will exhibit little distortion. Such a tube characteristic is pos- 
sessed, for example, by the 6L6 and 807 beam power tubes. 

As in the case of the previous constructions for the single-ended tube, various de- 
grees of complication can arise. For example, if the grids are driven positive so that 
grid current flows, the signal voltage at the grids will be distorted, and this distortion 
must be determined separately before the above construction can be concluded. An- 
other case is that where the mid-branch plate supply has an internal resistance that 
is adequately by-passed for the even harmonics, all except the d.c. component. This 
represents a particularly difficult problem that can be solved only by a series of ap- 
proximations. 

6. Reactive circuits. The previous circuits contained only resistances, linear or 
nonlinear. If reactances were present, such as the center-tapped output inductance, 
they were assumed infinite in value and so situated in the circuit as not to have any 
appreciable a.c. components flowing in them. However, many nonlinear circuits con- 
tain reactances of finite value that influence the behavior of the circuit directly, and 
hence must be taken directly into account. 

Owing to lack of space, only the case of an inductance in series with a nonlinear 
resistance and an a.c. source will be discussed here. Consider the circuit shown in 
Fig. 10. Here a source of a.c. voltage e is in series with a nonlinear resistance r and 
inductance L. The voltage e is a known function of time, and the terminal character- 
istic for r and the value of L is given. It is desired to find the current flow in this cir- 


cuit. 





1 See, for example, A. Preisman, Graphical constructions for vacuum tube circuits, McGraw-Hill Pub- 
lishing Co., New York, 1943. 
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We have the fundamental relation 
; di f 
et) =~ ar+L—- (7) 

dt 


Expressing Eq. (7) in terms of finite increments, we obtain 
; ‘ Ati 
é: + Aey = (4; + Aty)r +L a . (8) 
t 


In Eq. (8), it is assumed that at the start, the voltage e has a certain value e, the 


current 7 has a certain value 7:, and t=0. These are the initial conditions. During a 
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small time interval Af, e: is assumed to change instantly to e;+Ae, and remain at this 
value during the interval] At, and similarly i; is assumed to change instantly by an 
amount Az; and to remain at the value i,;+Az, during the time At. This is of course an 
approximation, sufficiently close if At is taken sufficiently small. Under these condi- 
tions Eq. (8) holds. 

The quantity L/At has a finite value if At is finite. It can represent the cotangent 
of some angle 6. Then—as far as Az; is concerned—the circuit consists of two resist- 
ances in series: that of r at the value 7,, and that of L/At. The graphical construction 
then takes the form shown in Fig. 11, where OA represents the initial value e,, and 
AB the initial current 7:. We now suppose that the voltage changes from e; to e,+Ae 
in a small chosen time interval] At, and let OD represent e,+Ae, so that AD repre- 
sents Ae;. 

The voltage across LZ is due to the change of current Az; and not due to 4, itself, 
which has already been established in L. This is indicated by the fact that OA=e, 
represents the drop across the nonlinear resistance 7 ; there is no voltage drop across L 
for 7; at the time ¢=0. Hence, in view of the above, a point C is located in line with B 
and directly over D, and through C line EC is drawn to represent L/At such that 


cot X ECB = L/At. (9) 


The line EC has been designated by the author as a finite operator because it re- 
sembles the Heaviside operator Lp. The intersection of this finite operator with the 
terminal characteristic of r in E gives the value of Az, namely, EJ. Here BJ repre- 
sents the additional voltage drop across r (in addition to the original voltage drop OA 
owing to 7,), and JC represents the voltage drop across L. In short, O0A+ BJ represents 
(i; +Ai;)r; JC represents L(Ai;/At); and OA+BC therefore represents e:+Ae,, and 


hence satisfies Eq. (8). 
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The point E is projected over to F directly above C and D, and FD represents 
then the new value of current 7;+Az;, at the end of the time interval At. Another 
small time interval can now be chosen, preferably equal to-the previous one, so that 
L/At remains at the same angle to the e-axis as before. We suppose that in this new 
time interval, e changes from e,+Ae, to e,+Ae,+Ae.. Letting OG represent the new 
value of voltage, we project F over to H directly above G. Through H we draw HK 
parallel to CE, intersecting the terminal characteristic for ry in K. Then KL represents 
the new increment of current Ai., EL the additional voltage drop across r, and LH 
the new voltage drop across L. It is evident that Eq. (8) is once again satisfied. It is 
also evident that IG represents (7; +A7,;+Az2), the new value of current at the end of 
the second time interval. 

Points B, F, and I represent three points on the overall terminal characteristic for 
L and r in series for the given function e(#). If e(t) is a periodic voltage, the overall 
terminal characteristic will spiral around counter-clockwise and ultimately form a 
closed curve, the steady-state solution for the given circuit and given function e(f). 
The initial open branches of this spiral represent the transient solution. If r is a 
linear resistance so that its terminal characteristic is a straight line instead of the 
curve shown in Fig. 11, the closed loop will be an ellipse inclined to both axes; if on 
the other hand 7 is nonlinear, the closed loop will be some form of distorted ellipse 
depending upon the nonlinearity of r. It can be shown from the graphical construction 
that the tangents to the closed loop at the points where it intersects the terminal char- 
acteristic for r are parallel to the e axis and hence perpendicular to the 7 axis. 








os eae. 








Fic. 12. 


7. Relaxation oscillator. Similar methods can be developed for r in series with a 
condenser C, and for LCr circuits, and for parallel as well as series arrangements. 
Owing to lack of space these will not be treated here.? An interesting case is that of a 
nonlinear resistance having a suitable negative branch, in series with a pure induct- 
ance. For graphical purposes the simplest form for the terminal characteristic of r 
is possibly that of three intersecting straight lines, as shown.in Fig. 12. Such a char- 
acteristic may be approximated by a tube having positive feedback, by a dynatron, 
etc. Usually a d.c. polarizing voltage is required, but this merely represents a transla- 
tion of the axes and does not materially change the construction or results as obtained 
in Fig. 12, in which the impressed voltage is assumed to be zero. 


2 Cf. Preisman, loc. cit., p. 109. 
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Suppose that the initial conditions at t=0 are that e=0, and that 1=AB, the 
peak current for the left-hand portion of r. Then C will be the starting point, where 
CO=AB. Through C the finite operator L/At is drawn corresponding to a time in- 
terval At. If AZ is sufficiently small, L/A? will be practically a horizontal line through C. 
In Fig. 12, L/At has been drawn with a finite tilt to clarify the construction, and is 
represented by CD. This finite operator curve intersects the terminal characteristic 
for r in D, as shown. 

The current therefore decreases from CO to DG. Point D is projected over to the 
i axis as point E. From E, EF is drawn parallel to CD under the assumption that the 
second time interval is equal to the first. The current now decreases from DG to FH. 
Point F can now be projected over to the 7 axis and the process repeated. It is clear 
from the figure that the intersections will proceed down the right-hand branch of r 
to I, hop over from I to J, directly opposite I, then proceed from J up to A, hop over 
to D, and repeat the first set of intersections. As At approaches zero, the finite opera- 
tor curve approaches a horizontal position, DG=CO=AB, and the points of inter- 
section become more and more closely spaced so that they form essentially all the 
points of ID and JA. 

The overall terminal characteristic is by definition all the points between C and 
K in that the overall impressed voltage has been assumed zero, so that the points 
must lie along the 7 axis, and the current range is from C to K. However, a more sig- 
nificant terminal characteristic in this case is the relationship between the current 
and the voltage across either circuit element. The voltage across the inductance, for 
example, is equal and opposite to that across r when taken in a circuital direction, 
since the algebraic sum of the two must equal the impressed voltage, which is zero. 

According to this definition, the terminal characteristic is represented by such 
points as D, F, etc.; in this case, it is lines DI, IJ, JA, and AD, traversed in the order 
given. This means that for the circuit given, the terminal characteristic is very simply 
given by a quadrilateral involving the two positive resistance portions of the terminal 
characteristic for r contained between their peak values A and I. 

The time required to traverse these portions depends upon the relaxation time 
for L in series with the incremental resistance of r for each portion, under the proper 
initial conditions. The time required to traverse the horizontal portions AD and IJ 
is infinitesimal, and is independent of the shape of the negative resistance portion AI 
provided it has no maxima or minima exceeding or less than A and I, respectively. 
The device operates continuously as an oscillator with a period of oscillation deter- 
mined by the two relaxation times. 

Similar conclusions can be drawn for. shapes of 7 other than three straight lines. 
For example, r can have the form of a cubic parabola. This case has been treated 
analytically by Van der Pol.* However, he started with an LCr parallel circuit or 
double-energy condition. For such a circuit the terminal characteristic is a closed 
curve or loop that exceeds the above quadrilateral in size. As C approaches zero, the 
loop shrinks and appears to have as its limit the above quadrilateral. However, the 
analytical method required that some capacity be present even in this limit, relaxa- 
tion case, and it has been suggested that in a practical circuit there would always be 
some residual stray capacitance present. 


3B. Van der Pol, The nonlinear theory of electric oscillations, Proc. I.R.E., 22, 1051-1086 (1934). 
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There are other graphical methods for handling the double-energy case, notably 
that by Liénard‘ and another by Kirschstein.® Unfortunately, these constructions be- 
come indeterminate in nature as C approaches zero, so that although the relaxation 
condition is suggested by them, it cannot be conclusively shown to be the limit form. 

The construction given here starts out with merely L and r, and requires no C 
for its argument. It appears to give the limit case directly and presents no indetermi- 
nate considerations. It has seemed to the author that the necessity for requiring a 
capacity to be present, no matter how small, was an unnecessary restriction, and that 
the argument advanced that any practical circuit would have some capacity, ap- 
peared to be rather irrelevant, since the notion of a circuit is in itself an idealization 
of what is really a field problem. In treating an electrical problem as a circuit problem 
one assumes that the circuit elements are ideal inductances or capacitances or resist- 
ances and develops the various theorems on this basis. 

Similar results can be obtained for a capacitance in series with a nonlinear resist- 
ance having an S-shaped terminal characteristic provided that it is turned through a 
right angle from that shown in Fig. 12, i.e., provided that it is a single-valued function 
of the current rather than of the voltage. A familiar example is the neon tube relaxa- 
tion oscillator employed to generate a saw-tooth voltage. It is also possible to develop 
a graphical construction employing the finite operator method for an LCr circuit, and 
in this case L or C may be permitted to approach zero, depending upon the position 
of the S-shaped characteristic for r, without the construction becoming indeterminate. 
For example, the construction reduces to the form given in connection with Fig. 12 
if C is made to approach zero and r has the terminal characteristic shown in the figure. 

8. Conclusions. This concludes the discussion on some graphical methods for solv- 
ing nonlinear electrical circuits. Simple series circuits involving resistance elements 
only, are very simply solved by finding the intersections of their terminal characteris- 
tics. This can then be extended to more complicated resistances in which the current 
is a function of two voltages, as in the case of a triode tube. 

The next circuit considered is that of the ideal balanced amplifier having perfectly 
matched tubes and feeding the load resistance through an ideal transformer. Here the 
coupling of the two tubes through this ideal transformer requires a special construc- 
tion involving the sliding of a rule at a fixed angle along the tube characteristics. The 
wave shape of the output and of the mid-branch currents is then discussed, and it is 
shown that owing to the symmetry of the circuit the former can contain only odd 
harmonics; and the latter, even harmonics. 

Finally, a simple case of a reactive circuit involving a nonlinear resistance in series 
with an inductance is treated. Here the concept of a finite operator curve correspond- 
ing to L/At is developed and this curve is employed to solve the circuit. Similar meth- 
ods are available for capacitive circuits and for double-energy circuits involving both 
L and C. The method is applied to a suitable negative resistance in series with an 
inductance, and it is shown in a direct manner that this circuit can produce relaxa- 
tion oscillations. 





‘A. Liénard, Etude des oscillations entretenues, Rev. Gen. Elec. 23, 901-946 (1928). See also 
P. LeCorbeiller, The non-linear theory of the maintenance of oscillations, Journal IEE (London) 79, 361-378 
(1936). 

’ F. Kirschstein, Uber ein Verfahren zur graphischen Behandlung elektrischer Schwingungsvorgange, 
Arch, Elek. 24, 731 (1930). 
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PRESSURE FLOW OF A TURBULENT FLUID BETWEEN 
TWO INFINITE PARALLEL PLANES* 


BY 
P. Y. CHOU 
California Institute of Technology 


1, Introduction. The solution of the Navier-Stokes differential equations for the 
steady laminar flow through a channel or a circular pipe is well known for its mathe- 
matical simplicity. The reason for this simplicity is that for such flows Prandtl’s 
boundary layer equations hold rigorously for the entire region of the fluid. In other 
words the boundary layer extends up to the center of the channel, whereas in the 
case of the flow around a solid obstacle there is only a thin layer of viscous fluid at- 
tached to the surface of the obstacle. 

The steady turbulent flow through a channel or a circular pipe is more compli- 
cated in the sense that all the equations of mean motion and the equations of double 
and triple correlation previously developed!:? have to be utilized to account for the 
mean velocity distribution in the entire region of the channel, and that they can not 
be further simplified by physical arguments as proposed, for example, by the bound- 
ary layer theory. However, if we examine the algebraic equation that represents the 
mean velocity distribution across the channel, we notice that it has functional be- 
haviour similar to that of the formula for the mean velocity distribution within a 
turbulent boundary layer.’ In other words, the turbulent flow in a channel bears 
some resemblance to the corresponding laminar flow on the whole, though its detailed 
structure is much more complicated as will be seen soon. 

In what follows we shall first determine the mean velocity distribution based upon 
the equation of mean motion and the equations of double correlation, by giving the 
triple correlations their values in the middle of the channel. This procedure leads to 
good results in the theory of the spread of turbulent jets and wakes (references at 
the end of II), but in the present case it only agrees with the experiment in the central 
portion of the channel, while it fails when the side is approached. We shall also see 
that the mean squares of the three components of the velocity fluctuation agree quali- 
tatively with observation in the corresponding region. 

The second determination given below for the mean velocity distribution utilizes 
equations of mean motion and both the equations of double and triple correlation by 
neglecting terms involving quadruple correlations. It will be shown that the triple 
correlations which represent the transport of turbulent energy play a particularly 
important role in the vicinity of the wall of the channel, and therefore can not be 
dispensed with for a better representation of the mean velocity distribution. 

From this second determination we shall find that neglect of terms involving quad- 
ruple correlations is justifiable as a first approximation. In other words the equations 





* Received Dec. 19, 1944. 

1 P. Y. Chou, Chin. Journ. of Phys. 4, 1-33 (1940). This paper will be referred to hereafter as I. 

2 P. Y. Chou, On velocity correlations and the solutions of the equations of turbulent fluctuation, Quart. 
of Appl. Math. 3, 38-54 (1945). This paper will be referred to as II. 

*N. Hu, The turbulent flow along a semi-infinite plate (unpublished). 
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of mean motion and of the double and triple correlations are sufficient in treating 
turbulent flow problems even though there is a wall present. Hence the mathematical 
procedure is comparatively simple in another sense that the building of equations 
satisfied by higher order correlations can be dropped up to the present degree of ap- 
proximation. 

The first determination reveals that the variation of the mean squares of the tur- 
bulent fluctuation is slower than the corresponding variation of the mean velocity 
distribution across the channel, which agrees qualitatively with experiment. In view 
of the fact that measurements of the mean squares of the components of turbulent 
fluctuation have not been reported systematically in the literature for the flow under 
consideration, we shall omit the quantitative comparison of the theory with the ex- 
perimental data now available on these quantities. 

In the second determination the mean velocity distribution will be calculated by 
assuming constant mean squares of turbulent velocity components across the channel. 
This is justifiable due to the slower variation of these functions across the channel, 
and furthermore the mean velocity distribution remains practically unchanged in the 
major portion of the channel—with the exception of the immediate neighborhood of 
the wall—when the constant values assumed for these functions are different from 
each other. This procedure of assigning constant values to the mean squares of the 
velocity fluctuations and then calculating the mean velocity distribution can be con- 
sidered as the initial step in a method of iteration which will be explained in §3 below 
in greater detail. 

In the final section we shall indicate the uncertainties connected with the correla- 
tion integrals pointed out before (II, §8). They are probably not important for the 
mean velocity distribution, because they involve possibly the mean squares of the tur- 
bulent fluctuation which are taken to be constant for the present calculation These 
uncertainties could be removed, if we had better experimental information on the var- 
iation of the turbulent level across the channel and on the velocity correlation between 
two distinct points in flows such as the one examined here. In other words the present 
theory is perhaps sufficient so far as the mean velocity distribution is concerned, and 
it points out the possibilities for future investigations in turbulence along both ex- 
perimental and theoretical lines. 

2. Mean velocity distribution based upon the solution of the equations of mean 
motion and of double correlation. As before (I, §4) we take the positive x-axis (x =<") 
as the direction of mean motion of the fluid, the y-axis (y=x*) perpendicular to the 
two parallel planes forming the channel, and the z-axis (2 = x*) parallel to these planes. 
The plane in mid-channel is chosen as the xz-plane. From the equations of mean 
motion we have 


2 
r12/p = — Uo — vdU/dy, (2.1) 
where 
2 
— 0p/pdx = U,/d, o = y/d. (2.2) 
The quantity 2d represents the width of the channel and U, is the so-called friction 
velocity. 
Equation (2.1) defines the shearing stress 712. in terms of y and dU/dy. Except in 
the immediate neighborhood of the wall the viscous stress is small, so 712 is a linear 
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function of y. On the wall —vdU/dy should be equal to U?, and r12 should tend toward 
zero as a limit. 

The components 723 and Ts; vanish due to symmetry, as pointed out before (I, §4). 

From now on for simplicity we shall neglect the action of viscosity in the form of 
laminar stress in all the equations of motion. A physical condition mentioned previ- 
ously is that all average values over time are functions of y only. Furthermore in 
the present special case in accordance with the definitions in II, Eqs. (5.3), those 
components of the slowly varying tensors a>,;, and b;, which have a single appearance 
of the index 3 either in 7 or k must be both identically zero. The vanishing of these 
functions is based upon the same argument as in the case of 723. The non-vanishing 
equations of the second order correlation (II, (8.2)) then become 











ii a eons au 4. @ wen, = — Gi aU — 5, + Rad (k — 5)q? — avk ant (2.3) 
p dy dy! a. a . 1 ie a eaaics 
1 dU , dU 2vk 
eee -— WW; = — do12 — — br. — —— W1W2, (2.4) 
p dy dy dy ? 
Oe andi ig ene Ee. O89 
dy * «ay ~ ga? zs 2 
d — dU 2 ye 2vk _ 
dy WeW, = — 2133 ae ae bss +- 32 (k —- 5)q’ =~ 2 ws (2.6) 


These are obtained by giving i and k the sets of values (1, 1), (1, 2), (2, 2), (3, 3), 
respectively. 
In the above four equations g is the root-mean-square of velocity fluctuation de- 


fined by 





q’ = w;w', (2.7) 
and k is a constant. The slowly varying tensors dni. and b;, should obey the diver- 
gence relations [II, (5.4) | 

2111 + G2122 + A2i33 = O, bir + bee + bss = 0. (2.8) 


The equation of vorticity decay [II, (7.11)] satisfied by Taylor’s scale of micro- 

turbulence A becomes in the present case 
— 14Gw,wdU/dy — 70Fq?/3\/3 = — 2vEq?/3X*, (2.9) 
where E, F and G are regarded as constants. 

The constant G in (2.9) is probably not important, for in the center of the channel 
the term involving G in (2.9) is zero due to the vanishing of dU/dy there, and in the 
immediate neighborhood of the wall w,w, vanishes although —dU/dy is large. Hence 
for simplicity we choose G to be zero. In fact the presence of G would only change our 
results slightly, as will be seen. The physical meaning of neglecting G in (2.9) is that 
the term that represents the creation of vorticity by deformation of the mean mo- 
tion is negligible when compared with those due to transport and decay. 

If G were set equal to zero, Eq. (2.9) yields 


\g/v = V3E/35F = Ro, (2.10) 


where Rp is a constant number, 
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Now we shall substitute the values of the triple correlations at the center of the 
channel into Eqs. (2.3)—(2.6) according to their odd and even properties as functions 
of y. We can set 


oo 8 _— 3 aa 3 —-— 3 
ww, = U,a\o, ws = Ua, wow? = U,a30, wi ws = Um, (2.11) 


where a, @, @3 and ay are four dimensionless constants. Although the factor U? in 
(2.11) is introduced for dimensional reasons, it is possible that these four constants 
are all independent of the Reynolds number of the mean flow. 

If we substitute from (2.10) and (2.11) into (2.3), (2.5) and (2.6), add the three 
together and take into account the conservation relations (2.8), we find that the mean 
square of the velocity fluctuation gq? satisfies the relation 





. & 2e dU 
— = «- = =], (2.12) 
u* 10R, U, de 
where 
a=a,;+aret+asz, R, = U,d/v; (2.13) 


R, is called the friction Reynolds number. 

Relation (2.12) is very significant, for it tells us that for large values of —dU/do, 
g? varies as the square root of —odU/de. Within a large portion of the channel, 
dU/de is proportional to a, so the dependence of g? upon @ is fairly linear. This linear 
dependence has been observed to some extent for w? by Wattendorf and Kuethe* 
and by Wattendorf and Baker,’ and has been anticipated in the light of von Karman’s 
law of similarity. 

If G were different from zero, the above procedure would lead to 


v =_R | 2(1 + 35G) — “| TOR ( os ~)" (2.14) 
— = c— 35G) —- — /10R,( ¢ — — — 
-_ U. de /s = 


T 





which has a functional behaviour similar to that of (2.12) for large values of —dU/de. 

It is apparent that Eqs. (2.3), (2.5) and (2.6) will determine w:, w2 and ws; sepa- 
rately. Here we encounter the uncertainty pointed out in II, §8 that the slowly vary- 
ing functions Gpmix and b;, may contain powers or even more complicated functions 
of g as factors, and the existing experimental data do not provide enough evidence 
for a quantitative comparison with these theoretical formulae. If, for the sake of 
mathematical convenience, we assume }j;, bo. and 633 to be constant, and de111, @2122 and 
d233, Which are odd functions of ¢, to be proportional to ¢, then w?, w2 and w? will 
behave very much like g?, that is, when @ is near zero, w?, w? and w? are constants, 
and when ¢@ is large and near unity, they are all proportional to the square root of 
—oadU/de. 

The equation for determining the mean velocity distribution is given by (2.4) 
which can now be written, on account of (2.10) and the condition that w,w, is con- 
stant, in the form 


‘ F. L. Wattendorf and A. M. Kuethe, Physics 5, 153-164 (1934). 
5 Th. von Kérm4n, Proceedings of the Fifth International Congress for Applied Mechanics (Cam- 
bridge, Mass. 1938), p. 349. 
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A ft 1 dU 2kR,a? 
—— i + a2) —— ——- = - peas (2. 35) 
U? , do Ri0} 


here we have set the odd function },. equal to zero for simplicity. 

As pointed out before the dependence of d2112 upon g in the above equation is also 
not known. If w? and g? are both regarded as constants, the mean velocity distribu- 
tion according to (2.15) is parabolic, which agrees with experimental data fairly well 
for the range of o from 0 to 0.8, and fails near the walls of the channel. This parabolic 
law of velocity distribution has been suggested by Stanton® in his measurements of 
flows through a circular pipe of which the channel is a limiting case. 

It has been calculated, though details will not be shown here, that this parabolic 
distribution of the mean velocity for constant g? and w, is not essentially changed if 
we solve for w?, w2, w2 and dU/do simultaneously under the further assumption 
that both ajq and b,;, are equal to zero. This condition is equivalent to the vanishing 
of (@,;w%+@,4w;:)/p, which means that the shearing interaction between the pressure 
gradient and velocity fluctuations is zero; it has been used in jets and wakes, as men- 
tioned before. The reason why the velocity distribution is parabolic even for this 
more rigorous treatment is not difficult to see without going into detailed calculations. 
For in the neighborhood of ¢=0, both w2 and g? are constants, so dU/de is propor- 
tional to ¢. When the values of o are near unity, both g? and.w? are proportional to 
the square root of —odU/de, and hence mutually proportional; consequently Eq. 
(2.15) again shows that dU/da is proportional to o even in the vicinity of the channel 
wall. We should anticipate, by the same argument, that similar simultaneous solu- 


tions for w?, w2, w2 and dU/de would hold true even under the more general condition 
that bu, be and 633; be constants and G11, G2122 and d233 be proportional to o as men- 








tioned previously. 
In §4 below we shall compare the numerical values of Ro, R, and a of (2.12) with 


available measurements. 
3. Equations of triple correlation and the mean velocity distribution. The non- 


vanishing equations of the triple correlation [II, (8.3) ] for the present problem can be 


written in the form 




















dU or dU 3 dry 
3u* we. — — wwe = — be1111 — — C111 + = ea (3. 1) 
dy dy ' dy p” dy 
inatiianin dU 4 d —— b dU 4 1 (2 dr 12 4 <=) 3 2) 
20,08 — + — wwe = — Danis — — Crs + — | 2712 — T ; (3. 
wW1W, dy dy 12 111 dy il 2? 1 dy 11 dy 
ee BU ——_ dU 1 dry. dr22 
ws —-+— ww oe bori22 — C122 + 7, i eae + 2712 ’ (3.3) 
dy dy . dy p” dy dy 
: dU dU 1 dr12 
wow? — + — wiw2u? = — be1133 —— — C133 + ee ’ (3.4) 
S dy dy . dy p? dy 
i —< dU 3 dT» 
— we = — beie%2 —— — Coxe + — T22 —> (3.5) 
dy dy p? dy 


* T. E. Stanton, Proc. Roy. Soc. London. (A) 85, 366-376 (1911). 








1945] PRESSURE FLOW BETWEEN PARALLEL PLANES 203 





d —_os dU 1 dro 
dy — =" boie3s en — Co33 + > T33 rs , (3.6) 


These are obtained by giving i, k and / the sets of values (1, 1, 1), (1, 1, 2), (1, 2, 2), 
(1, 3, 3), (2, 2, 2), (2, 3, 3), respectively. The other component tensor equations in 
which the index 3 appears an odd number of times, namely, (1, 1, 3), (1, 2, 3), (2, 2, 3) 
and (3, 3, 3), are all identically zero, as are the corresponding equations of the second 
order correlation. 

From the discussions in the previous section it is apparent that Eqs. (2.3), (2.5) 
and (2.6) are used to determine the mean squares of the fluctuation components, and 
elimination of the triple correlations between these three equations and (3.1), (3.3) 
and (3.4) respectively will give a more accurate determination of them. As pointed 
out before, existing experimental data are not accurate enough to give a quantitative 
comparison. with the theory, and we shall not go into these detailed calculations here. 
Furthermore Eqs. (3.5) and (3.6) lead to quantities which are still beyond experi- 
mental proof; discussions of them will also be omitted for the present. 

The elimination of the triple correlation ww? between (2.4) and (3.2) leads to the 
equation for the mean velocity distribution. Before writing down this equation we 
shall introduce a few more simplifications. In the first place the even function be, 
which may depend upon g as mentioned previously (II, §8), is assumed to be a con- 
stant; likewise the odd function ¢,,2 is taken to be proportional to y and is put in the 





form, 
ens = 2cU.e/d. (3.7) 


It is also possible that the dimensionless number c may be a function of g and there- 
fore an implicit function of the coordinate y. 

The quadruple correlation w?w? in (3.2) is of the same order of magnitude as 
(w;we)? and w? w2. As a first approximation we shall neglect all of these terms and it 
will be shown afterwards in §4 that this approximation is justifiable. In short, (3.2) 
defines the triple correlation ww: approximately by 

a 4 dU 
ww? = — dban2e — cU,o / —-: (3.8) 
do 

Utilizing the above relation and (2.10) which is derived from the equation of 
vorticity decay, we find, after setting by. in (2.4) equal to zero for mathematical con- 


venience, that 





a dU d U,o 
as aarmes aig Siem Serene ae, (3.9) 
U, do do dU 
da 
where 
a= —(w+tan)/Ur, 6 = 2kR,q /RoU;. (3. 10) 


The physical meaning of the three terms in the above equation is as follows: the term 
in a represents the creation of turbulent energy partly due to deformation of the mean 
flow (I, §3(a)) and on account of d21.2. partly contributed by the shear due to the pres- 
sure fluctuation (@ w+ .w;)/p; the term in } denotes the decay of turbulence; the 
term in c denotes the transport of turbulent energy. 
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The definitions of a and b in (3.10) show that they depend upon g? and w?, and are 
therefore functions of a. Since g? as well as w? as shown in the previous section varies 
much more slowly than dU/do itself across the channel, we shall regard them as 
constants as the initial step to solve for dU/de. This initial process can also be re- 
garded as the first step in the method of iteration in solving the present problem of 
turbulent flow. The second step will be to substitute this expression obtained for the 
mean velocity into (2.3), (2.5) and (2.6) after eliminating the triple correlations by 
means of (3.1), (3.3) and (3.4), and to solve for w?, w? and w?. As the third step in this 
procedure, we utilize these values of the mean squares of the fluctuation components 
and solve (3.9) again for dU/do, and see whether the new result agrees with the solu- 
tion obtained in the first step. Obviously this procedure of obtaining alternately the 
mean velocity and mean squares of the turbulent fluctuation can be extended in- 
definitely. 

In the present paper we shall not follow this refined method of approach; instead 
we shall solve (3.9) by assigning constant values to a, b and ¢, or to g? and w?, and com- 
pare the different solutions by varying these constants. The result will be that except 
in the immediate neighborhood of the wall of the channel, the different mean velocity 
distributions according to (3.9) for the different sets of a, b and ¢ respectively agree 
well with each other and with experiment, showing that the variation of the mean 
squares of the turbulent fluctuation across the channel does not influence the mean 
velocity distribution very much. 

The solution of (3.9) with constant a, b and c is 


do? = aU/bU, + cAseUleUr + Ao, (3.11) 
where A; and A? are two constants of integration. 

If A; in (3.11) is zero, then (3.10) gives a parabolic law of velocity distribution 
and a must be negative, since } according to its definition in (3.10) is positive. The 
presence of the term in A, gives the so-called “logarithmic law” of velocity distribu- 
tion which holds true especially in the neighborhood of the wall of the channel. Hence 
the product cA; can not be zero. It is apparent that this exponential term in U/U, is 
due to the presence of the triple correlation in Eq. (3.9). 

The boundary conditions used to determine the constants cA,, Az and the ratio 
a/b are: 


when oc =0, U=U.-: when o = 1, U=0, —dU/U,do = ~. (3.12) 


The value U, denotes the maximum velocity of the flow in mid-channel. We have 
chosen the derivative —dU/U,do on the wall of the channel to be infinite. In fact, 
it should be R, which is a fairly large number. Since we are interested in the mean 
velocity distribution within the channel proper, substituting infinity for the friction 
Reynolds number R, gives a good approximation. 

The boundary conditions (3.12) render (3.11) into the following final form, 


(e*VelUr — KU./U, — 1)o® = — «(Ue — U)/U, + eel s[1 — eee VIUr], (3.13) 


where x=)b/c and b/a=2[exp (xU./U,) —xU./U,—1]/k. 

Equation (3.13) expresses the mean velocity defect (U.— U)/U, as a function of o 
with two parameters x and U,/U,. The presence of these two constants may appear 
at the first sight to contradict the experimental velocity defect law formulated by 
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von Karman,’ according to which (U,— U)/U, should be independent of the Reynolds 
number of the mean flow which, in turn, is a function of the ratio U./U,. A close ex- 
amination of the experimental data shows, however, that this discrepancy is not 
serious. In the first place von Kaérm4n’s velocity defect law can only hold true in the 
central portion of the channel and there is a dependence of the velocity defect upon 
the Reynolds number in the vicinity of the channel wall. It has been shown that for 
flows in circular pipes U,/ U, increases from about 19 to 33 when the friction Reynolds 
number 2aU,/v changes from +/105 to 10%, 2a being the diameter of the pipe.* In the 
second place even formula (3.13), which does indicate the dependence of (U.— U)/U, 
upon U,/U,, can only account for the mean velocity distribution in the interior of 
the channel for a given set of constants a, 6 and c in (3.11), and these constants have 
to take another set of values in the turbulent boundary layer on the wall, although 
the same functional behaviour of (3.11) still prevails within the layer.’ This point will 
be discussed in greater detail in the following section. 

The quantity U./U, in (3.13) is given by experiment; then the constant x is fixed, 
for instance, by passing the theoretical curve through the experimental point at 
o =0.7. In view of the variation of the ratio U./U, with the Reynolds number of the 
mean flow, we shall choose a few different values of x and calculate the mean velocity 


TABLE 1. (U.—U)/U,; 























(1) (2) (3) (4) (5) 
: et | Obs. —0.1 0.0 +0.1 0.2151 

0.00 0.00 0.00 0.00 0.00 0.00 
0.10 0.16 0.06 0.05 0.05 0.05 
0.20 0.38 0.23 0.22 0.21 0.19 
0.30 0.66 0.52 0.50 0.48 0.44 
0.40 1.10 0.95 0.91 0.88 0.81 
0.50 1.64 1.50 1.47 1.42 1.34 
0.60 2.33 2.22 2.19 2.15 2.08 
0.70 3.13 3.13 3.13 3.13 3.13 
0.80 4.28 4.31 4.38 4.50 4.75 
0.90 6.30 5.92 6.17 6.60 7.72 
0.93 — 6.58 6.91 7.51 9.30 
0.96 8.81 7.39 7.88 8.74 11.84 
0.98 — 8.12 8.76 9.93 15.01 
0.99 — 8.64 9.40 10.82 18.21 
1.00 - 9.86 | 10.94 13.07 2 





distribution. This will lead to different values of U./U,. But we shall see that for all 
these cases the mean velocity distributions agree with each other and with experiment 
within the channel proper. 

Let us calculate the mean velocity distribution for the values of x equal to —0.1, 0, 
0.1 and 0.2151, and determine the corresponding values of U./U, by passing the theo- 
retical curves through the experimental point at ¢=0.7. The equations that deter- 
mine (U,—U)/U, for x= —0.1 and 0.1 are given respectively by 





7Th. von K4rmén, Proceedings of the Fourth International Congress for Applied Mechanics, 


Cambridge 1934, p. 70. 
8S. Goldstein, Modern developments in fluid dynamics, vol. 2, The Clarendon Press, Oxford, 1938, 


p. 338. 
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o? = 0.2785(U- — U)/U, — 1.039[e9-1Ue-U/Ur — 4], (3.14) 

o? = — 1.5870(U. — U)/U, + 2.662[1 — e-9-1(Ue-¥) 10s], (3.15) 

For the case x=0, we can get a limiting equation by letting x approach zero in 
(3.13), or we can solve (3.9) directly by setting 6 equal to zero. The latter procedure 


leads to 


(U. — U)/U, = 10.94[1 — (1 — o)!/2], (3.16) 


tl 


where 10.94 is the value of (c/a)*/?. 
The case when x=0.2151 is represented by 


2 


c= 1 ‘ais e@79-2151(Uc—-U) /U ge. (3.17) 


This is the solution of (3.11) with a set equal to zero; the numerical value 0.2151 


stands for the ratio b/c. 
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Fic. 1. Velocity distributions in a channel. 
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The experimental values of (U.—U)/U, which are taken from a paper by Gold- 
stein® are given in column (1) of Table 1; the corresponding theoretical values ac- 
cording to (3.14), (3.16), (3.15) and (3.17) are tabulated in columns (2), (3), (4) and 
(5), respectively. 

From this table we see that as the value of x increases from —0.1 to +0.2151, 
U./U, changes from 9.86 to ». Hence 0.2151 is the maximum limiting value of x 
for constant a, 6, c in (3.9). Equation (3.17) shows in this limiting case that the mean 
velocity distribution in the whole channel is “logarithmic.” 

In order to avoid confusion, only the solution (3.14) for x= —0.1 is plotted in 
Fig. 1. The circles represent Dénch’s measurements!® found for U,d/v equal to 
8.710‘, U,, being the average value of U over a cross section of the channel. The 
crosses reproduce Nikuradse’s results" for U,d/v equal to 3.3104. It is seen that 
apart from the immediate neighborhood of the channel wall, agreement between 
theory and experiment is satisfactory. 

4. Relation between the present theory and some known experimental data. From 
the foregoing calculations we see that we can subject to experimental test not only 
the mean velocity defect distribution (U.—U)/U, and the mean squares of the 
fluctuation components, but also the relation (2.10) between \ and gq and the relation 
(3.8) which approximates the triple correlation w,w2. Let us examine relations (2.10) 
and (2.12) first. 

The experimental data used by Taylor in his statistical theory’ are, in c.g.s. 
units; U.=114 cm/sec, U,=5.39 cm/sec, p=0.00123, »v=0.14, d=12.3 cm. Since 
according to (2.10) Ro=Aq/v is a constant, we can compute Ry from the values of 
\ and gq in the center of the channel. In Taylor’s table \? in mid-channel is equal to 
2.9 cm’, so \ is 1.7 cm. The mean magnitude of the velocity fluctuation g at this 
point is roughly 1.2U, [cf. (4.2) below]. These values then give 


Ro = 78.5, R, = 474. (4.1) 





We have shown above in §2 that w? behaves very much like g?. The experimental 
values of w?/U? across the channel, determined by Wattendorf and Baker® for the 
flow with Reynolds number 109,000, can be represented by 


wi/U; = 0.412(1 + 27.20)". (4.2) 
As far as the order of magnitude is concerned, w? can be put equal to g?/3. By com- 
paring (4.2) with (2.12), we find that 
RoVa/V10R, ~ 3 X 0.4 = 1.2. (4.3) 
If we use a parabolic representation of the mean velocity distribution that goes 
through the experimental point at ¢=0.7, we have 
— 20dU/U,do = 25.60". (4.4) 
Then @ in (2.10) becomes 25.6/27.2 =0.94. Putting this value of a in (4.3), we obtain 
R?/R,~10X1.4/0.94=15, while from (4.1) we find that R3/R,~13. 


® S. Goldstein, Proc. Roy. Soc. London (A) 159, 473-496 (1937). 

10 F, Dénch, Forschungsarbeiten des Ver. Deutsch. Ing. no. 282 (1926). 

it J. Nikuradse, Forschungsarbeiten des Ver. Deutsch. Ing. no. 289 (1929). 
122 G. I. Taylor, Proc. Roy. Soc. London (A) 151, 456 (1935). 
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This shows the order of agreement between the two sets of values obtained from 
two entirely different experimental sources. It must be pointed out, however, that 
the number R, =474 for the experimental value of \ used in the above calculation 
may be too low for the flow to be in a fully developed turbulent state. 

Equation (2.10) also shows the dependence of Ry upon the quantities E and F 
which occur in the definitions of the double and triple correlation functions between 
two distinct points [II, (6.8), (6.11) ]. The measurement of these functions separately 
will give another check on the value of the number Rp discussed above. 

We next study the values of the three constants a, b and c in Eqs. (3.9) and their 
physical significance. 


(1) x= —0.1. According to (3.10) 6 must be positive, so c in this case must be 
negative. The definition of c/a from (3.13) and (3.14) gives 
c/a = 20/0.2785 = 72. (4.5) 


Hence a must be also negative. If dey2 in (3.10) were zero, a becomes of the order 
of 0.4 and c is equal to 29. From the definition of c in (3.7), we find that 


tan » 50U,/d, (4.6) 


which is 29 times larger than 27,2dr12/p*dy, a term of the same order of magnitude as 
the one in the quadruple correlation dw?w2/dy. Hence all these terms are negligible as 
a first approximation. 

After the value of ¢ is known, the triple correlation function w,w? is determined 
uniquely according to (3.8); the constant term 4be11. is fixed by the value of w,w? at 
o=0. 

By means of b= —0.1c, the definition of 6 in (3.10) and the values of R,/R% and 
q'/ U? given before, we find the numerical value of k to be of the order of 16. 


(2) x=0, (c/a)'/?=10.94, c/a=120. This gives results similar to those in case (1) 
and consequently the terms in the quadruple correlations in (3.2) are still negligible. 
In this case c and a must be negative as in the foregoing example. The meaning of 
b=0 is that the term due to the decay of turbulence in (3.9) is negligible when com- 
pared with the other two. 


(3) x=0.1. In this case ¢ should be positive. Then the definition of c/a from (3.13) 
and (3.15) gives c/a=20/1.5870~13, and a must be also positive. The condition 
that a be greater than zero changes the picture a great deal, for then @2112 in (3.10) is 
negative and its magnitude is greater than w?. Nevertheless the terms in the quadruple 
correlations are still negligible, if the absolute value of de2 is, say, a few multiples 
of w?. 

(4) «=0.2151, a=0. Here we have b/c =0.2151. Putting this value into the defini- 
tion of 6 from (3.10), we have c~10XkX1.2/15, which is about 12, if & is of the 
order of 15. This gives cy. of (3.7) equal to 24U%a/d, a quantity still about 10 times 
greater than the terms involving the quadruple correlations in (3.2). The physical 
significance of a=0 means that in (3.9) the term due to deformation is small when 
compared with the terms due to transport of and the decay of the turbulent energies. 
Obviously from the equations of double correlation (2.3)—(2.6) the magnitude of k 
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can be determined by measurement of the mean squares of the velocity fluctuation 
components. 

From the four alternative cases discussed above we see that although the measure- 
ment of the mean velocity distribution alone will not single out which one is the 
correct theoretical mean velocity distribution, measurements of the variations of the 
higher order correlation functions across the channel will decide this question. For 
example, experiment on the triple correlation w,w? in (3.8) will decide whether c is 
negative or positive, and the theoretical pattern for the mean velocity distribution 
can thus be determined. 

From another angle the above four special cases can also be considered to repre- 
sent the mean velocity distribution in four parts of the channel. In the central 
portion, we have negative a and negative c [cf. (3.10), (3.8)]. Since both ¢ and dou 
can be functions of the coordinate a, their values may change at the various points 
of the channel. It is possible that a may eventually become positive as o increases 
near the wall, while c which was negative in mid-channel, increases to zero and finally 
becomes positive on the wall of the channel. 

According to its definition in (3.10), b is positive and is a monotonically increasing 
function of the distance from the center; likewise x increases monotonically with ¢, 
if ¢ has already become positive. This increasing property of x as the wall is ap- 
proached is substantiated experimentally. In Hu’s theory of the turbulent flow along 
a semi-infinite plate,’ the mean velocity distribution in the turbulent boundary layer 
can be represented by an equation analogous to (3.11), and the value of b/c is equal 
to 0.4 instead of 0.1 as in case (3). Hence our present solution for mean motion only 
covers the channel proper; if the boundary layer on the channel wall is approached, 
the solution should be replaced by Hu’s result. In fact it is well-known experimentally 
that the turbulent boundary layer on the wall covers the region 30 <R,(1—a) <250. 
In Dénch’s measurement’® cited above, R, is equal to 3630, so the range 0.931 <¢ 
<0.992 represents approximately the turbulent boundary layer on the wall and we 
should expect formula (3.13) to fail in this region. A rigorous theory to explain the 
mean velocity distribution for the entire channel including the boundary layer 
might not be impossible according to present indications, but the actual mathematical 
manipulation involved would be much more complicated than that in the present 
treatment. 

5. Conclusion. Based upon the foregoing analysis in the cases of the four values 
of x for the motion of a turbulent fluid through a channel, we may conclude that the 
velocity defect distribution (U.—U)/U,, which is practically independent of the 
Reynolds number of the mean flow within the channel proper according to von 
K4rmAn, is also independent of the magnitudes of the turbulent fluctuation when the 
flow has reached the steady turbulent state. The question as to whether the above 
conclusion can be generalized to state that the double and triple correlation distri- 
butions across the channel when expressed in terms of the frictional velocity U,, 
namely, the ratios w,w,;/U? and ww,w,/U?, are also independent of the Reynolds 
number of the mean flow and of the correlations of still higher orders remains to be 
seen theoretically as well as experimentally. In any event, the friction velocity U, 
probably plays an important role for turbulent flow problems invoiving the presence 
of a wall, as in the present problem. 
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DIFFUSION IN TURBULENT FLOW BETWEEN 
PARALLEL PLANES* 


BY 


J. C. JAEGER 
University of Tasmania 


1. Introduction. The equation of diffusion in a turbulent fluid 


0*x 1—2p 0x ox 
—~4—~ 2 (1) 


dz” Z 02 «Ox 


~ ~ 





(x stands for temperature, vapour concentration, or whatever property is being stud- 
ied, x is measured in the direction of mean flow and 2 in the perpendicular direction, 
and »p is a constant determined by the degree of turbulence of the fluid) was intro- 
duced by O. G. Sutton! and extensively studied by W. G. L. Sutton,? who considered 
a number of cases of diffusion in the semi-infinite region z>0. It has been shown by 
Pasquill* that for the semi-infinite region the theory is in good agreement with ex- 
periments, both on evaporation and on heat transfer. 

In this note a number of results for symmetrical flow in the finite region 0<z<2l 
will be given; it is assumed that 2/ is small enough for the power law velocity profile 
to hold up to the centre of the region. Such cases are of some practical interest, and 
may provide an indication of the behaviour to be expected in the much more difficult 
problem of heat transfer in a circular pipe. Also they are interesting generalizations 
of known solutions of the equation of conduction of heat in the rod 0<z<2l, with 
constant temperature, or flow of heat, at its ends. 

The method used will be that of the Laplace transformation. W. G. L. Sutton 
(loc. cit.) remarks that if p=1/2, equation (1) reduces to the equation of linear flow 
of heat, and he gives a treatment of (1) which is a generalization of Goursat’s treat- 
ment of the equation of conduction of heat. It is well known that the Laplace trans- 
formation method is particularly well suited to the solution of specific problems in 
conduction of heat, and that its advantage increases as the complexity of the problem 
increases. This suggests that the method may have the same advantages when applied 
to (1), and, in fact, this proves to be the case. All the results of W. G. L. Sutton’s 
paper can be obtained more shortly in this way, and explicit expressions for the solu- 
tions for more complicated boundary conditions, composite regions, etc., can also be 
derived. 

In Section 2 the standard problem of evaporation in the semi-infinite region is 
solved as an illustration of the method, and for comparison with later results. In 


* Received Dec. 26, 1944. 
10. G. Sutton, Wind structure and evaporation in a turbulent atmosphere, Proc. Roy. Soc. (A), 146, 


701 (1934). 
2 W. G. L. Sutton, On the equation of diffusion in a turbulent medium, Proc. Roy. Soc. (A), 182, 48 


(1943). The notation used here is that of this paper, except that the symbol E£ is introduced in (4). The 
variables x and z in (1) are dimensionless quantities defined in Sutton’s paper. 
3 F, Pasquill, Evaporation from a plane, free-liquid surface into a turbulent air stream, Proc. Roy. Soc. 


(A), 182, 75 (1943) 
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Section 3 two other results for the semi-infinite region are given for completeness. In 
Sections 4-6 the most interesting cases of symmetrical flow in the region 0<z<2/ 
are studied. The solutions given here are formal only, but in all cases they may be 
made rigorous by the verification process described elsewhere.‘ 

Equation (1) has to be solved in the region x >0, and in a prescribed region of z, 
with boundary conditions in x 


x—>x(z), as x—> + 0, (2) 
x finite, as x—> ©, (3) 


In all the problems considered below x (z) will be zero, that is the temperature or 
vapour concentration in the fluid is zero in the plane x =0. 
There are also boundary conditions in z, which will be expressed either in terms 
of x, or of 
ax 
E = — Bz'-*? — . (4) 
dz 
This quantity E is the local rate of diffusion across the plane z=const., and B is a 
known constant (defined by Sutton, Joc. cit.) involving the fluid and its degree of 
turbulence. 
The constant p in (1) is restricted in Sutton’s theory by the inequality 0<p<1/3, 
and we assume here 0<p<1. 
With the substitution 


x = 272, (S) 
(1) becomes 
an 1 a2 Pp ao 
—+— — -—-—2-— =90. (6) 
02? 2 OZ 2? Ox 


Introducing the Laplace transform of 2 with respect to x, namely 


ce) 


o* = f e"Qdz, (7) 
0 
we obtain from (6) the subsidiary equation for 0*, 
a oe AS “ 
, dz? Zz dz 2? 


2. The semi-infinite region z>0. Boundary conditions: x (z)=0, z>0. x=xo, 
constant,’ when z=0, x>0. x finite, as sz «, x >0. Here (8) becomes 


ao* 1 do p? 
——+— ——-—(s+—)m®=0, <«>0, (9) 


. 
dz? z dz z 


4H. S. Carslaw and J. C. Jaeger, Operational methods in applied mathematics, Oxford, 1941, §58, 
and J. C. Jaeger, Radial heat flow in a circular cylinder with a general boundary condition, Proc. Roy. Soc. 
N.S.W., 75, 130-139 (1942). 

5 For the procedure see, e.g., Carslaw and Jaeger, loc. cit. 

* For shortness, boundary conditions will usually be written in this way; it is implied, of course, that 
x—xo as z—> +0 for fixed x>0. 








212 J. C. JAEGER [Vol. III, No. 3 


to be solved with 
x* = z?90* > xo/s, as 2+ 0, (10) 


and 
(11) 


x* finite as 27 ©, 


The solution of (9) which satisfies (11) is K,(zs!), and since 





2°K,(zst) > 2>"T'(p)s-#” as z—>+0, (12) 
it follows that the solution of (9), (10), and (11) is 
s*2*? 
x* = xo —— s!'K,(zs!). (13) 
I'(p) 
Now it is known that? 
si?-1K ,(zs!) is the Laplace transform of evar f e~“u?—'du. (14) 
22/42 
Thus the required solution is 
1 = mare 
Xx = x0 e~“uP—"du. (15) 
['(p) J 242 
This is the result given by Pasquill (Joc. cit., (9)). 
It follows from (15), or directly from the transform E* of E, that 
Bxo2!-??x-? 
i; i 4 > Oo (16) 
I'(p) 
and 
z Bxo2'—??x'-? : 
f Edx — ——————_» as 2—+ 0. (17) 
0 (1 — p)T(p) 


3. Two other results for the semi-infinite region z>0. The results to be derived 
here are both for the case x (z) =0, and x(z) finite as z>@. 
If the boundary condition at z=0 is: E—+Eo, constant, as z—>+0, the solution is 





Eoz?? a 
x = _ — f e~“u-?—"du. (18) 
2BI(1 — p) 22/42 
Also 
Ey2??-1x? 
as z—>+0. (19) 





een >? a 
Bor (1 — Pp) 
This is proved exactly as in Section 2, using (14). 
7S. Goldstein, Operational representations of Whittaker's confluent hypergeometric function and Weber's 
parabolic cylinder function, Proc. London Math. Soc. (2) 34, 104 (1932), (15) and (24). Alternatively the 
result can be obtained by the use of the inversion theorem for the Laplace transformation, subsequently 
deforming the line integral into the contour (— ©, 0+); cf. Carslaw and Jaeger, loc. cit., §39. The same 


remark applies to the derivation of (18) and (21) below. 
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If the boundary condition at z=0 is 


Ox 
hx — z!-2P eg hxo, (20) 


2 


where h and x» are constants, the solution is 





X02” renee J (uz) + au??J_,(uz 
9 ef raw get = a pus) du, (21) 
2?"T(p) Jo 1 + 2au?” cos pr + a*ut? 


where a=I'(1—p)2'-**/hT (p). 

To prove (21) the inversion theorem for the Laplace transformation, (24) below, 
has to be used, and the line integral must be deformed into (— «, 0+). 

The result (15) was derived for a constant value of x on the boundary z=0. The 
solution for the case in which x is a prescribed function of x on z=0 can be obtained 
from (15) by Duhamel’s theorem. The same remark applies to the cases of Sections 4 
and 6. Correspondingly, the solutions of the problems of Sections 3, 5 with E a pre- 
scribed function of x on z=0 can be obtained in the same way. 

4. The region 0<z<l. x(z) =0. x =xo, constant, when z=0, x>0. E=0, when 
z=l, x>0. This corresponds to the region 0<z<2/ with flow symmetrical about 
z=l, and with x =xo on z=0 and z=21, for x>0. Thus, for example, it gives the solu- 
tion of the problem of heat transfer from the parallel planes z=0 and z=21, both 
maintained at constant temperature xo, and with symmetrical flow between them. 

Here we have to solve (9) with boundary conditions (10) and 


E*=0, when 2z=l. (22) 
By (12) the solution of (9) which satisfies (10) is 
xoz?2!~? 
x* = —_— s!1K, (2s!) + Az°I,(zs4). 

I'(p) 
The unknown A is found by substituting in (22), and we have finally 
x02?2'-Pstr—1[ K ,(z54)T p_1(Is*) + Ip(2s!) K p_(Is) | 

T'(p)I p-a(/s}) 


x is found from (23) by using the inversion theorem for the Laplace transforma- 
tion [cf. Carslaw and Jaeger, loc. cit. | 





x* = (23) 








1 ytio 
x =— e**x*(s)ds (24) 
2rid +i 
ef — e**s47—1[K (284) I p_1(Is4) + Ip(2s!) Kp-1(Is') Jas (25) 
~ Oeil (p) T,-x(is!) 


y¥—t0o 


where y>0. 
The integrand of (25) is a single valued function® of s. It has a simple pole at s=0 


of residue 
2>-13-7T (p), (26) 


®G. N. Watson, Theory of Bessel functions, Cambridge University Press, Cambridge 1922, §3.71, 
(17) and (18) 
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and simple poles at s= —a?/I?, where +a,, r=1, 2, - - - , are the zeros (all real and 

simple®) of 

J p-1(@) = 0. (27) 
It is easy to show [cf. Jaeger, loc. cit.] that the line integral in (25) is equal to 27i 

times the sum of the residues at the poles of its integrand. Evaluating these we get 


finally 7 
X02? 2 atte ol! J ,(sa,/l) 


- ———— (28) 
2>-*P°T(p) m1 J*(c,) 








x = 7 


The most interesting quantity is the value of E as z-0. Either from (28), or di- 
rectly by calculation of its transform, this is found to be 








BxoV (1 aime ) ce) en~arz/l* 7 a Gr) 
E- = P > i-o as z—-+ 0. (29) 
I'(p)2??-2]?? nos ai—2PJ (ar) 
Also, as z—>-+0, 
z Bx [2-2p 2 
f Edx + ———_ 9, (x/I'), (30) 
0 2(1 — p) 
where at 
2 T'(2 — p)2%-2? 2 emarzll Jy_ (ar) 
$, (s/f) =1- p = . (31) 


'(p) r=1 at 2PJ (ar) 


For small values of x//? the value of pe Edx given by (30) reduces to the value (17) 


for the semi-infinite region, and 


ae 23-27 ¢/§*)!-» 
Se —————- (32) 
r'(p) 


In Fig. 1 graphs of these quantities are shown for p=1/9, the value commonly 
found in wind tunnel experiments. Curve I shows the result (32) for the semi-infinite 
region, and Curve II the value of ¢9)(x//?) given by (31) for values of x//? for which 
the difference between (31) and (32) is important. For larger values of x//? than those 
shown the exponentials in (31) are negligible. 

In the case of heat transfer the quantity (30) gives the amount of heat taken up 
from the region 0 to x of one of the planes. 

5. The region 0<z<i/. x(z)=0. E— Eo, constant, as z—-0, x>0. E=0, z=1, 
x >0. This corresponds to the region 0<z<2/ with flow symmetrical about z=/, and 
with constant diffusion across the planes z=0 and z=21. 

Here, proceeding as in Section 4, we find 


: Eoz2? { Ip_1(Is*)K p(zs4) + Kp—s(Is#)Ip(zs') } (33) 
io BI(1 — p)st?+47,_,(Is!) 





x 


The most interesting quantity to evaluate in this case is the value of x as z—+0. 


This is found to be 
2E,(1 — pl?” x 
B Ee + =| (34) 

*G. N. Watson, loc. cit., §§15.25, 15.21. For the method of calculating their values in practice see 
John R. Airey, Bessel functions of small fractional order and their application to problems of elastic stability, 
Phil. Mag. (6), 41, 200 (1921). 
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where — 22”-11(p) os e287 1 (cep) 
§{3/) ° (35) 
2p(2— p) T(2— p) ni a??tJ_,(a,) 
and the a, are the roots of 
Ji~p(a) = 0. (36) 


For small values of x//?, (34) tends to the result (19). 
If p=1/9 it is found that the difference between (34) and (19) is less than 1% 
for values of x//? up to 0.3, while for greater values of x//* the exponentials in (35) are 
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almost negligible. In the case of heat transfer, (34) gives the surface temperature of 
one of a pair of planes to which heat is supplied at a constant rate per unit time per 
unit area, and which are cooled by turbulent fluid flowing between them. 

6. Two cases of symmetrical flow in the region 0<z<2/. First let us consider 
boundary conditions x (z) =0, and 


X = xo, constant, when z=0, x >0 (37) 
E=0, when z= 2l, «> 0. (38) 


Here the regions 0 <<z</ and /<z<2/ must be treated separately. We write xi(21) 
and £;(z;) for the values of x and E in] <z<2l as functions of 2; = 2/—z in this region. 
The boundary conditions at the surface of separation z=z,=/ are 


X= x, (39) = — &;, (40) 








216 J. C. JAEGER [Vol. III, No. 3 


A solution of (9) which satisfies (37) is 


X02) PstP-! 
x* EEE 2K ,(zs}) oa Az?I,(zs}), 
I'(p) 


and a solution of (9) with z replaced by 2; which satisfies (38) is 


* _ 
x1 = C2PI_,(z:54). 
The unknowns A and C are found by substituting in the transforms of (39) and 
(40), which gives 
* 


m x02'?apst?-47_,(z,5!) 
'— P(p)[T-p (ts!) pa(Is) + Tp(ls!)Tip(Is4) | 


(41) 








with a rather longer formula for x*. As in Section 4, x and x: are evaluated by the use 
of the inversion theorem and the results are 


xo2!-P2iT'(1 — p) e~F47!"gP-lT (cep) J p(atr) J —p(zar,/l) 














x1 = X0— (42) 
la Ip ~ J3-1(cr) + F3 (ar) 
Xo2!-?g? = e~207/U gyp—27 (za,/l) 
[= a= ~—— e : 2 (43) 
1°T (p) r=1 J2_1(a;) + J? (a,) 
where the a,, (r=1, 2, - - - ), are the positive roots!® of 
J_p(a)J p1(a) — Jp(a)Ji_p(a) = 0. (44) 


. . . . . . . z= ~ rr . . 
As in Section 4 the most interesting quantity is the value of J Edx as z—0. This is 


found to be 








fz Bxol?-*” (2), ’) (45) 
tas —> —————— z/i), 5 
x11 - Pp op (%, 
where 
2 2*-*7(1 — 9) © e-eerit 
o> (x/l) = 2 -—-—_—_—__—— >} —_—_ (46) 
‘ [T(p)]? 7 at-*”{ J? (a,) + J2_,(ar) } 


For small values of (x//*), (46) behaves like (32). Its value for p=1/9 and for val- 
ues of (x//?) for which the difference from (32) is important is shown in Curve III 
of Fig. 1; for larger values of x//? the exponentials in (46) are negligible. 

The result (45) gives the evaporation from the region 0 to x of the plane z=0 if 
there is no flow over the plane z= 21. 

Finally we consider the case in which the boundary conditions are x‘ (z) =0, and 


X = x0, constant, s=0, x>0, (47) 


and 
x =0, when z= 2l1, x>0. (48) 
Here, proceeding as before and writing x; for the value of x in] <z<21/, E, for the 
value of Z in this region, and z;=2]—z, we find 
10 In the case k=1/2 when (1) becomes the equation of linear flow of heat, (44) becomes cos 2a =Qand 
similarly (51) and (52) become sin 2a=0. 
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x= x0(1 - =") a a it (49) 
2») r(p) SL at ?I2_,(a) Bt-7J2(B,) 
x1 = Axo (=)" xo "a (ae é es » (50) 
l VI(p) a1 a2-?J?_, (az) 62-? J? (B;) 
where the a, are the positive roots of 
J (a) = 0, (51) 
and the 8, are the positive roots of 
J p-1(B) = 0. (52) 


This problem is that of heat transfer between the plane z = 0 at constant tempera- 
ture xo, and the plane z=2/ at zero temperature, by turbulent fluid flowing between. 
The quantity of heat taken up from the region 0 to x of the plane z=0 is determined 





by 
* Byxol??? 2 
lim f Bha  eaemie tall), (53) 
Powe we 2(1 — p) 
where 
(3), 4.2 2p(.1— p)x 1+ p+ 2p 
op (x/l) = ——— a 
§ 2(1 + p) 
23-2p 1 ane C) enzarl i? e~78r/07 
[T(p)]? tar bat*9J2_,(ar) BS??? (B,) 


and the a, and 8, are defined by (51) and (52). For small values of (x/FP), ¢® (x/P) 
behaves like (32). For larger values its behaviour for the case p=1/9 is shown in 
Fig. 1, Curve IV, and for still larger values the exponentials in (54) are negligible. 

The quantity of heat taken up by the region 0 to x of the plane z=2/ is deter- 
mined by 


Bxol??? ay 














“3 2 
— ; = / 
jim, E,dx 1 op (x/l), (55) 
where 

wk 2A p)x (1+ 2~)(1 — p) 
6, (z/1) = ————__- 

P 2(1 + p) 

23-27(1 — p) x { enzarll? en Belt? | ae 

[r(p)]? Fat bet*7J2_ (az) BE*”J2(B,)) © 


A portion of the curve of ¢%(x/P) for p=1/9 is shown in Fig. 1, Curve V; for 
larger values of x//? the exponentials in (56) are negligible, and for smaller values than 
those shown in the figure o continues to decrease rapidly. 








218 


ON THE STABILITY OF TWO-DIMENSIONAL PARALLEL FLOWS 
PART II.—STABILITY IN AN INVISCID FLUID* 
BY 


C. C. LIN 


Guggenheim Laboratory, California Institute of Technology 


7. General considerations. The criteria of Rayleigh and Tollmien. At the end of 
Part I, we have shown that the study of the stability problem in an inviscid fluid 
gives valuable information provided it is kept in mind that we are actually dealing 
with the limiting case where the Reynolds number becomes indefinitely large. The 
study of the stability of two-dimensional parallel flows in an inviscid fluid is usually 
regarded as being quite complete, through the work of Rayleigh and Tollmien. Their 
results show that instability depends very much upon the occurrence of a point of 
inflection in the velocity profile. However, it seems that physical interpretations of 
such general results are not well known. Such an interpretation will be given in §§9, 10 
of this part. There are also several points in the mathematical theory which require 
further development and clarification. These will be brought out for further considera- 
tion in §§7, 8. 

We now proceed to make a critical survey of some aspects of the stability problem 
in an inviscid fluid. First, let us summarize the conclusions obtained by Rayleigh and 
Tollmien. These can be conveniently described as the necessary and the sufficient 
conditions for the existence of a disturbance, self-excited, neutral, or damped. 

1) Necessary conditions for the existence of a disturbance. 

a) If the flow possesses a self-excited or neutral mode of disturbance with finite 
wave length, the velocity profile has a flex at some point y=y,, where yi<y.<yo. 
Furthermore, in the case of a neutral disturbance, the phase velocity must be c=w(y,). 

b) If the flow possesses a damped mode of disturbance, no immediate conclusion 
can be drawn. 

2) Sufficient conditions for the existence of a disturbance. So far, the sufficient condi- 
tions are known only for symmetrical and for boundary-layer velocity distributions. 
The results may be stated as follows. 

a) There is always the neutral disturbance given by c=0, a=0, $(y) =w(y). 

b) If w’’(y,) =0, for y1$<y.<ye, there is a neutral disturbance with c=w(y,); fur- 
thermore, if w’’’(y,) #0, self-excited disturbances also exist. 

Discussion. The condition w’’’(y,) #0 involved in (2) (b) will be shown to be ac- 
tually unnecessary, by an improved method of proof to be discussed in the next sec- 
tion. The statement in (1) (6) regarding damped disturbances differs from the origi- 
nal conclusion of Rayleigh and Tollmien. Indeed, in the work of Lord Rayleigh, the 
solution is taken to be valid all along the real axis. Hence, in accordance with the dis- 
cussion of §5, Part I, such considerations do not include damped disturbances. How- 
ever, Rayleigh and Tollmien did not distinguish between an amplified disturbance 
and a damped disturbance, because they regarded them as complex conjugates. As 
pointed out in §5, this is not permissible. In fact, if we accept the original conclusions 


* Received May 18, 1945. Part I of this paper appeared in this QUARTERLY, 3, 117-142 (1945). 
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of Rayleigh and Tollmien, a profile without a flex could not execute any kind of dis- 
turbance. This can hardly be reconciled with our intuition regarding the state of 
affairs in a real fluid at infinitely large Reynolds numbers.' According to the present 
interpretation, only damped solutions can exist. Such a conclusion is also borne 
out by the investigations for a viscous fluid.? It is to be noted that the neutral and the 
self-excited disturbances, existing under the condition w’’(y,)=0, are free from the 
effect of viscosity inside the fluid, because the neutral solution is also regular at y=y, 
where w=c. Hence, we may conclude that disturbances essentially free from the eff ect 
of viscosity inside a fluid can exist only for velocity distributions with a flex. 

The results of Rayleigh and of Tollmien discussed above tend to give the impres- 
sion that the occurrence of a flex in the profile is the decisive factor in the determina- 
tion of instability not only in the case of an inviscid fluid, but also in the case of a 
viscous fluid.* However, the investigation in Part III will show that this is by no 
means the case. When instability first occurs, as one increases the Reynolds number, 
viscous forces still play a dominant role, and the main characteristics of the behavior 
of the fluid with respect to a disturbance do not depend upon the occurrence of a flex 
in the velocity curve. Indeed, it is physically improbable that a slight change of the 
pressure gradient in the case of a boundary layer—which may cause a change from a 
velocity curve without a flex to one with a flex—should cause a radical change in the 
essential characteristics of stability. As we shall see later, the instability of a boundary 
layer depends more on the outside free stream than on the occurrence of a point of 
inflection. It might be argued that the free stream is analogous to a point of inflection 
in that a vanishing curvature is involved; but even if this is admitted, we must still 
note that the essential features in this case are not obtained from an investigation 
neglecting the effect of viscosity. Indeed, from inviscid investigations, it is concluded 
that a boundary layer with zero or favorable pressure gradient is stable, except for 
the very trivial type of disturbance with infinite wave-length and zero phase velocity. 
The present investigation shows that all boundary-layer profiles can be unstable, and 
exhibits results in agreement with the physical suggestion just discussed. 

It thus seems that any conclusion obtained from inviscid investigations must not be 
taken over directly to the case of the real fluid, where the stability phenomenon is largely 
controlled by the effect of viscosity and not decided primarily by the occurrence of a flex 
in the velocity curve. 

Indeed, even when we are mainly interested in the behavior in the limiting case 
of infinite Reynolds numbers, the existence of a flex is not as significant as it may 
appear to be at first sight. The existence of neutral or amplified disturbances has so 
far been proved only for symmetrical and boundary-layer types of velocity profiles. 
This may not be true for other types of velocity profiles, e.g., when the walls are in 
relative motion. The following example will bring out this point. Let us consider the 
velocity distribution w(y)=A+Bsiny, yi<y<y2, which has a flex at y=0 if 
yi1<0<ys. According to the above necessary conditions, the only possible neutral 





1 This is the objection of Friedrichs, loc. cit. (Ref. [5]) p. 209. (The references are listed at the end of 
Part I.) It must also be noted that the non-linear terms are not negligible in the case of an ideal fluid. 
We shall consistently restrict the magnitude of our disturbances so that the effect of viscosity is always 
more important than the effect of non-linearity. 

2 See figures in Ref. [27]. 

? See Taylor's discussion on p. 308 of Ref. [70]. 
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disturbance is the one with c=A. Then the equation of disturbance (6.21) reduces 
to ¢’’+(1—a?)¢=0. It has the solution 


¢(y) =Csin {V1 — a? (y— y)}, 


which vanishes at y=¥4;. If @(y2) is also required to vanish, we must have 


V1 — a? (yo — 91) = ne, (n = integer), 
and hence 


w=1— [nm / (V2 -_ yi) |. 


Thus, if ye—3yi<7, there is no possible neutral disturbance; if ye—y1=7, there is the 
one with a=0; if t<ye—y,<2z, there is one with a0; in general, if mr<y2o—y1 
<(m-+1)z, there are m neutral disturbances with a0. In the last case, there are 
also m points of inflection in the velocity profile. 

It can thus be seen that the general shape of the velocity profile plays a very im- 
portant role even in the limit of infinite Reynolds numbers. Indeed, it will become 
clear from Part III that the eigen solution with eigen values a=c =0 is not as trivial 
as it might appear at first sight, for it actually represents a limiting case with R->o. 
This solution exists for symmetrical and boundary-layer profiles, but its existence is 
not immediately evident for other types of profiles. 

In spite of all these points against the decisive nature of the flex, it must be ad- 
mitted that its occurrence certainly makes the motion comparatively unstable. This 
can be expected from the original results of Rayleigh and Tollmien, and can be seen 
more clearly from the interpretation of the mechanism of inertia forces to be given 
in §§9, 10. However, these results must not be taken to indicate any decisive nature 
of a flex. The essential features of instability can only be obtained through considera- 
tion of the effect of viscosity. 

We shall now conclude this section by making some critical discussions of Heisen- 
berg’s classification of velocity profiles and the use of broken linear profiles for the 
study of stability problems. 

Heisenberg’s classification of velocity distributions. Heisenberg attempted the case 
of flow between solid walls in relative motion with the condition that Re(w—c) van- 
ishes only once for yi<y <2 (loc. cit., p. 592). Regarding a? as small, he approximated 
the condition (6.18) by Ki(c) =0 [cf. (6.26), (6.24) ]. He then classified the profile into 
four classes: (7) those for which Ki(c)=0 has a complex root; (41) those for which 
Ki(c) =0 has a real root; (iz) those for which the real part of Ki(c) vanishes for a cer- 
tain real value of c; (iv) those for which none of the above three cases is true. Heisen- 
berg concluded that the first class is unstable, the second generally unstable, the rest 
stable. 

In discussing the validity of these conclusions, the following point must be borne 
in mind. If we can show that a certain type of disturbance exists for a?=0 and 
aR— «, it may also be expected to exist for sufficiently large values of aR and suff- 
ciently small values of a?. However, the non-existence of a certain type of disturbance 
for a?=0 and aR—~ does not exclude the possibility of its existence for finite values 
of a? and aR. It appears therefore that we can only expect to conclude the instability 
of a velocity distribution by discussing the roots of Ki(c)=0. Thus, apart from some 
flaws in Heisenberg’s mathematical deductions, only the first two classes can have 


any decisive significance. 


1945] STABILITY OF PARALLEL FLOWS 221 


If Ki(c) hasa root witha positive imaginary part, the motion is unstable. If K,(c) 
has a real root, Heisenberg shows that the motion would be unstable when the effect 
of viscosity is considered. This will be studied more fully in a generalized form in §11. 
However, if Ki(c) has a root with a negative imaginary part, we cannot conclude the 
instability of the flow by taking the complex conjugate of K,(c)=0 (as Heinsenberg 





did). For if rg 
y, > Y, 
Jaw —c)?=0, % Yo 

then (cf. Fig. 5) 
f dy(w — c)~-? = — 2miRo, < 

c’ Fic. 5. Path around the critical point in the case c;<0. 
where R> is the residue of (w—c)-* at yo. In fact Ro= —w¢’ /we*. Now K,(é) is the 
complex conjugate of [edy(w—c)-*. Hence, 

K,(é) = 2riRy =a 2riws' /we?, 


which does not vanish unless w/’ =0. Hence, the equation Ki(c) =0 tells us nothing 
about the existence of the root ¢ or any other root with a positive imaginary part. 

Thus, Heisenberg’s attempt appears to be not as successful as Tollmien’s later 
work [75], which at least brings out the characteristic properties of symmetrical and 
boundary-layer distributions. A complete classification of velocity distributions, how- 
ever, is not yet existent. 

Approximation using broken linear profiles. Some investigations of Lord Rayleigh 
were carried out by approximating the velocity profile with straight-line segments. 
With this approximation, the solutions of (6.21) can be expressed in terms of elemen- 
tary functions. Lord Rayleigh also tried to verify his conclusions by considering the 
roots of Ki(c)=0, using the same approximation for the velocity. However, the re- 
sults of his investigations are doubtful, because the number of roots obtained for 
Ki(c) is equal to the number of corners chosen in the approximation. This was demon- 
strated by Heisenberg to be inherent in the method of approximation. The general 
idea is as follows. As discussed above, the stability condition (6.18) may be approxi- 
mated by Ki(c) =0 in certain cases. Although Rayleigh’s approximation may be made 
very close so far as the velocity distribution is concerned, the approximation to 
(w—c)~* is always bad in the neighborhood of the corners. Consequently, the integral 
K,(c) is not properly approximated. In fact, a continuous broken profile w(y) does not 
allow itself to be continued analytically to the complex y-plane without introducing 
discontinuities (cuts). It thus appears that all results deduced from the consideration 
of broken profiles must be regarded with reserve. The same criticism applies to 
Tietjen’s work with the viscous fluid. His analysis failed to give a minimum Rey- 
nolds number below which all small disturbances are damped out. 

8. Rigorous proof and extension of Tollmien’s result for the existence of unstable 
modes of oscillation. In this section, we want to give a rigorous proof of the existence 
of amplified solutions of (6.21) satisfying the second and the third boundary condi- 
tions of (6.22) when the velocity profile w(y) has a flex at y=y,, i.e., 


w''(y.) = 0. (8.1) 
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The idea of the proof is essentially the same as that used by Tollmien, but the method 

is improved. It has the further advantage of enabling us to extend the results to cover 

cases where w’’’(y,) =0,—a condition which had to be excluded by Tollmien. 
According to previous results, the neutral disturbance must have a phase velocity 


c equal to 

C. = Ww, = w(y,). (8.2) 
Let the corresponding value of a be denoted by a,. The essential idea of the proof is 
(1) to show that there exist eigen-values of c and a>0O in the neighborhood of the 
values of c, and a, such that the imaginary part of c does not vanish, and then (2) to 
show that the imaginary part is actually positive. The first statement can be expected 
and can be readily established, if we can show that the left-hand sides of (6.18)—(6.20) 
are analytic functions f(a, c) of the two variables a and c in the neighborhoods of a, 
and c,. For if this is true, we can always solve f(a, c) =0 for c as an analytic function 
of a, (there may be more than one branch), by the implicit function theorem. Hence, 
there is at least one value of c corresponding to every real value of a in the neighbor- 
hood of a=a,. Furthermore, by (8.2), this value of c, being unequal to c,, cannot be 
real, and the first part of our result is established. 

To prove the analyticity of f(a, c) seems to be a trivial problem. Nevertheless, we 
shall find below that it is impossible to establish it in the neighborhood of (a, c) = (0, 0). 
The chief problem in the proof is to overcome the difficulty caused by the singular 
point of the differential equation (6.21). 

If w—c +0, we can write (6.21) in the form 

id 
Ol Ghd ao cing me ©, (8.3) 
w—c 


Let us now consider a simply-connected region R of the y-plane which encloses the 


y-PLANE 


Ss c=W(y) 











Fic. 6. The region of analyticity of the inviscid solutions. 


points y=; and y=‘, but excludes the point y,, the passage from 4; to yz being taken 
in the lower half of the y-plane. We consider also a neighborhood S of y,, mutually 
exclusive with the region R. Let us regard the relation 


c= w(y) (8.4) 


as mapping the regions R and S into two regions R’ and S’ of the c-plane (Fig. 6). 
If the mapping is one-to-one, (as can be expected if w’(y)#0 for y1<y<y2), these 
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regions will also be mutually exclusive. Then, if we restrict y to R and c to S’, the 
coefficients of (8.3) are analytic functions of the independent variable y and the pa- 
rameters a and c. Hence, a fundamental system of solutions of (8.3), which we denote 
by o:(y; a, c) and ¢e(y; a, c), are analytic functions of the three variables y, a, and c. 
We understand that y is restricted to the region R, ¢ is restricted to the region S’, 
while a may be in any finite region enclosing a,. Thus, (for example), 

a _ a, 6) dely1; a, €) (8.5) 

oi (y2; a, ¢) 2 (ye; @, €) 

is an analytic function of the variables a and c, as we want to prove. 

We note that in the neighborhood of (a, c) = (0, 0), the above reasoning fails. The 
region R (which has to enclose the point y=) and the region S (which has to enclose 
the point where w=c=0) cannot be taken to be mutually exclusive. In fact, f(a, c) 
presumably has a singular point at the point a=0 (a logarithmic branch point). We 
shall discuss this case a little more closely at the end of this section. 

Let us proceed to show that there actually exist values of c=c(a*) with a positive 
imaginary part corresponding to positive real values of a. This is necessary because 
the usual argument of taking complex conjugates has been shown to be invalid. For 
this purpose, we consider the power series 


c= at(S)a-m+T(S)a-mwete (8.6) 
“as, “ 2thant/, : ; 
where \=a?.‘ Since ) is restricted to real values, the important point to be shown is 
that the first of the derivatives in (8.6) for which the imaginary part does not vanish is of 
odd order. Then, by taking values of X slightly greater or smaller than \,, we can al- 
ways make c;>0. For these values of c and a?, we can continue our solution ¢(y) 
analytically so that it is given along the real axis between 9; and ye, thus obtaining an 
inviscid solution. 
Let us now consider (8.3), writing \ for a?. We have 
w’’ 
Oe TEP > (8.7) 


Let ¢ be an eigen function with X, c as the corresponding eigen-values. Then 





Li ) pm ” 2 al {1 4 w’ = (8 8) 
“(pr) = or a*py ne iwce & dg, . 
where 
Op dd dc 
we? ie ae 


We distinguish two cases: (1) the point y=y, is a simple root of w’’(y) =0; (2) the 
point y=¥y, is a multiple root of w’’(y) =0. 
In the first case, w’’’(y,) #0. In the limit A-A,, c—c,, Eqs. (8.7) and (8.8) become® 


‘ Since d\/da #0 at a=a,, the correspondence between a and ) is one-to-one in the neighborhood of 
A, = a. 

5 A subscript s denotes that the parameters \ and ¢ are put equal to A, and c, respectively. A sub- 
script \ denotes differentiation with respect to A. 
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L.($s) = $8’ — Abs — ara ¢. = 0, (8.9) 
j_ Cs 


, w’’ w" de 
L (rs) == drs sid Ashrs ee Prs = {1 + (5 ) (8. 10) 
w Ge dx 


(w — c,)? 


From these, we deduce that 


ee ee Tie oie ie (2) he 
ds 8 Pre Prs 8 ds ad dy 1 DsPrs DrsPs om \ 2 ar io 


(w— 


Now, @) satisfies the same boundary conditions as @ does, because those conditions 
are satisfied by ¢ for each pair of values of \ and c, and ¢ is an analytic function of 
them. Hence, integrating ¢,L.(¢,s) —¢xsLs(¢,) between the limits (1, ye), we have 


dc v2 w’”’ 
(29) f* "ays fais 0 
A/a y, (w— C)? os 
iay/ f° wo, ‘(w — Cs *dy. 8.11 
(5)- af i oo y ( ) 


The denominator of the above expression is equal to 





or 


Vz—-Us 
f (wi! y + Fwiry? +--+ (why + dw!’ 9? +--+ 2G + Wuhay +--+ )dy 
Vvi—Vs 


wi'! : V2—-Us f 1 
alee if yy tAet A +--+ pay, 
w,* a 


Vi—-Veg 


where @,, is the value of ¢, at y=y,, and Ao, A1,: ~~ are real. Hence, the imaginary 
part of the above expression is 7¢,,w/’’ /w,*. Since ¢,(y) is real and ¢,, does not van- 
ish,* we have arrived at the required result. The above argument is a rigorous formu- 
lation of Tollmien’s work. 

In case w’’(y) has a multiple root at y=y,, the proof of Tollmien does not hold, 
but the above method can still be carried through. The restriction must be made, 
however, that the point y, is a point of inflection where w’’(y) actually changes its 
sign. Then, y, is a root of w’’(y) of odd multiplicity, and the first of the derivatives 
wi*(y,), w’(ys), -  - which does not vanish is of odd order. Such a point always exists 
when the curvature of the velocity curve has different signs at y; and ye. If we differ- 
entiate (8.7) m times with respect to A, we have the following equation for each value 


of n: 





 . d’ w”’ 

L(dxm) = ngyo—v + > C,dxn-»—( - (8.12) 
al a’ \w-e 

Let w’’(y) have the root y, up to the multiplicity 2m+1, m>0. Then, (8.10) is 

regular in a neighborhood of y=y,, and the value of (dc/dd), as given by (8.11) is 

real. Let us consider the boundary value problem of the differential equation (8.10), 

requiring $,, to satisfy the same boundary conditions as ¢,. The solution can be ob- 





 Tollmien [74], p. 92. 
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tained from @)(y) by making A—A,, and is moreover real along the real axis, by a 


direct consideration of (8.10). 
Continuing the same argument with equations of the type (8.12) with n=2, 3, 
- ++, 2m and c—c,, \—A,, we find that 


d% d3¢ d?™¢ 
(Gr )or (r™)ay + + + 5 (drome; (a), a. af Ga) 


are all real. Finally, for »=2m-+1, we obtain a relation of the type 


; dc 2m+1 v2 a ae q2mti¢ v2 . 
(2m + 1)! =) f w'o,(w — c) dy+ —- f ¢.dy = real. (8.13) 
8 vi Vu 








Just as in the case of the equation preceding (8.11), it can be easily seen that the above 
integral Siw!" (w—c)-em dy has the imaginary part w?"*®)g?, /(2m+1)! (wi )@mt® 
while the other term on the left-hand side of (8.13) is real. Thus, (d2"+'c/d)\2™*"), has 
a non-vanishing imaginary part. This is the result desired. 

This completes the proof of the existence of amplified solutions near the neutral 
solution ¢=c,, a=a, when the velocity curve has a point of inflection. 

The proof of the existence of amplified solutions near the neutral solution c=0, 
a=0 cannot be so easily formulated into a rigorous form. From the solutions (4.14), 
it is very easy to obtain the solution ¢; which approaches the eigen solution ¢ = w(y) 
as c—0, a?—0, with a? = O(c). The solution is 


uv 
or = — cwi w—o f (w — c)~*dy 
vi 


4 {1 + af ay(w - of ayw —cy?+-- , . (8.14) 


v1 v1 


As can be easily verified from (6.21) the condition that ¢,; be an eigen function is 


ve 
cw, + a f (w — c)ordy = 0. (8.15) 
vy 


1 


dc 1 ve 
(=) = —f wdy, (8.16) 
dx 0 wi vi 


and that the imaginary part of (d*c/dX?)> is 2r(dc/dd)ow{/w}’, which is positive if 


there is one flex in the velocity profile (w:’ >0). However, the real part of (d*c/dv*)o 
becomes logarithmically infinite, and hence the argument is not rigorous. Also, it does 
not seem easy to make suitable modifications and extensions in case w’’(y:)=0. It 
should be remarked that Tollmien’s proof is not essentially different from the argu- 
ment just given. 

Similar considerations can be applied to boundary-layer profiles, and similar re- 
sults can be obtained confirming and extending Tollmien’s original results. 

9. Physical interpretation of instability in an inviscid fluid. The fact that the in- 
stability of a two-dimensional parallel laminar flow is so closely connected with the 
occurrence of a point of inflection in the velocity profile demands a physical inter- 


From this, it follows that 
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pretation. Since Eq. (6.21) is essentially the vorticity equation, we would expect 
w’’=0 to indicate a maximum or minimum of the vorticity —w’ of the main flow. 
This is actually where the explanation is to be found. 

Since we have neglected the effect of viscosity, a fluid element maintains its vor- 
ticity throughout the motion. From this point of view, a two-dimensional parallel 
flow may be regarded as the motion of a large number of vortex filaments under the 
action of each other. Filaments of equal vorticity are arranged in the same layer, and 
the whole flow is built up of a collection of such layers. 

The following physical interpretation is based upon the fact that a fluid element is 
accelerated in such a field if it is associated with an excess or a defect of vorticity. 
These considerations were originally developed by von K4rm4n’ for the interpreta- 
tion of the failure of the simple vorticity-transfer theory of fully developed turbulence 
as applied to the case of parallel Couette flow. The idea is developed in greater mathe- 
matical detail here in this section and the next. It will be noticed that the considera- 
tion is essentially two-dimensional, and hence is even more suitable here than for 
fully developed turbulence, where the fluctuations are essentially three-dimensional. 
An alternative interpretation of the results of Rayleigh and Tollmien, but still based 
upon vorticity considerations, will also be given to demonstrate the role of the viscous 
forces. 

Let us imagine a disturbance of the flow such that an element £, of fluid of the 
layer L; is interchanged with an element £; of a neighboring layer Lz. For definiteness, 
let us suppose that the layer Z2 has a higher vorticity than the layer Z; in the undis- 
turbed state. Since E; preserves its vorticity, it will appear to have a defect of vor- 
ticity when it is in LZ». Similarly, Ez, appears to have an excess of vorticity. 

Let us fix our attention on one element, say E>. It will be shown in §10 that a fluid 
element with an excess of vorticity is accelerated in the direction of the positive 
y-axis with an acceleration I'-'/f}v’(x, y) }*¢0'dxdy, where {¢ (y) is the gradient of 
vorticity of the main flow, v’(x, y) is the component of the disturbing velocity per- 
pendicular to the direction of flow, and T is the total strength of the vortex filaments 
corresponding to the disturbance. Examining the signs of the various quantities in 
the acceleration formula, we can easily see that E, is accelerated toward a region of 
higher vorticity if the gradient of vorticity does not change sign anywhere in the fluid. 
Thus, E> is accelerated toward LZ». A similar consideration holds for the element £,. 
Hence, in either case, the fluid element is returned to the layer where it belonged (by the 
acceleration due to its interaction with other vortex filaments). The motion 1s there- 
fore stable when the gradient of the vorticity does not vanish. 

When there is an extremum of vorticity, an interchange of fluid elements on op- 
posite sides of the extremum does not give rise to an excess or a defect of vorticity. 
Furthermore, the gradient of vorticity vanishes there, and has opposite signs on 
opposite sides of that layer. It can easily be seen from the above acceleration formula 
that the restoring tendency mentioned above is largely impaired in such a case. Thus, 
exchanged fluid elements are not as strongly forced back by the action discussed 
above. Such an exchange constitutes a disturbance because there is an exchange of 
momentum. Thus, a disturbance may tend to persist and perhaps to augment. The 
motion is not necessarily stable. 





7 Cf. discussions of the vorticity transfer theory of turbulence in his general lecture at the Fourth 
International Congress for Applied Mechanics [19]. Some developments in that direction were continued 
by C. B. Millikan (unpublished). 


1945] STABILITY OF PARALLEL FLOWS 227 


The above discussion is based on very general considerations and does not depend 
on the consideration of a periodic wavy disturbance as used in the mathematical 
analysis. We shall now support the above argument by considering a neutral wavy 
disturbance, with the understanding that if such a disturbance can persist (except 
for the exceptional case of infinite wave-length and zero phase velocity), the motion 
is presumably unstable. From these considerations, the importance of viscosity in the 
inner friction layer will also be brought out. 

Let us consider an observer moving with the phase velocity of a neutral wavy 
disturbance. He will observe a stationary pattern of the flow (see Fig. 7).* Closed 
stream lines are inevitable unless the disturbance has no v-component of velocity in 
the critical layer w=c, for the 
flows on opposite sides of the criti- 
cal layer are in opposite directions 
relative to the observer. It appears at tt et 


unlikely that the v-component of 
the disturbance should be zero Ss <x > 
throughout that layer. Indeed, it 
has been shown to be impossible 
mathematically. Thus, whenever age aaa gg agg 


a neutral disturbance persists, it 








. ; << — 

involves a steady exchange of fluid 

elements on opposite sides of the 

critical layer. Fic. 7. Stream lines of a neutral disturbance as observed by 
If the effect of viscosity is to be an observer moving with the wave velocity. 


negligible, fluid elements on the 

same stream line must have the same vorticity. If the gradient of vorticity of the main 
flow is zero or small near the critical point, it is easy to compensate this small differ- 
ence of vorticity by the vorticity of the superposed flow, while the “scale” of disturb- 
ance [as measured in order of magnitude by u'/(du'/dy)| remain the same as that 
of the main flow. It is thus not impossible to find a neutral disturbance for which the 
effect of viscosity is negligible. The motion may be unstable. 

On the other hand, if the gradient of vorticity of the main flow is finite, the super- 
posed small disturbance must also give a finite gradient of vorticity. This means that 
the “scale” of the disturbance must be very small in the critical layer. The diffusion 
of vorticity by the effect of viscosity is then inevitable. It is thus impossible to find a 
neutral disturbance for which the effect of viscosity 1s negligible. The motion is inertially 
stable. 

10. Acceleration of vortices in a non-uniform field of vorticity. In the foregoing 
physical interpretation of inertial instability, we have considered the acceleration of 
an element of fluid in a two-dimensional parallel flow when this element of fluid does 
not have the same vorticity as the surrounding layer. We are now going to derive the 
explicit formula for the acceleration. The derivation shall be made in two different 
ways: (1) by kinematical considerations (using vorticity theorems) ; (2) by considera- 


8 This figure is due to Lord Kelvin (loc. cit.). He pointed out that the facts discussed here are “sur- 
prising,” but did not attempt to explain their connection with the mechanism of hydrodynamic stability, 

* This follows at once from Rayleigh’s original results, if we apply it to the region between this layer 
and the solid wall (cf. Tollmien, loc. cit., 1935). 








228 c.. &. La [Vol. III, No. 3 


tions of the pressure gradient. In either method, we shall consider a perfect fluid in 
accordance with the stability problem under consideration. 

1) First derivation, by kinematical considerations. For definiteness, let us consider a 
two-dimensional flow between two solid walls, which we shall take to be y= +b. Let 
the velocity components of the main flow be 


U = w(y), V = O, (10.1) 
and those of the secondary flow be 
u’ = u'(x, y), v = v(x, y). (10. 2) 
The distribution of vorticity of the main flow is 
fo = fo(y) = — w'(y), (10.3) 
and that of the secondary flow is 
dv’ du’ 
PP Se eee e's (10.4) 
Ox Oy 


The latter distribution shall approximate a vortex “at” the point (&, 70). Thus, if 
the signs of ¢’ and {> are the same (or opposite), we have essentially a small element 
of fluid having an excess (or a defect) of vorticity near the point (&, 70). 

The stream function for the secondary flow is 


1 
V(x4,9) = —- —{f o’(E, n)G(x, y; & n)dédn, (10.5) 
2 


with 
oy’ 1 ff (e 2) re) Gi b n)déd 
lille oe ’ —Ol\xX, ¥; &, ’ 
"teat ars y; &, n)dédn 


dy’ 1 P e 
v(x, y) = oe —ffz (&, 0) ht y; §, n)dédn. 


In these expressions, the integrals are extended over the whole region between the 
planes. The function G(x, y; —, 7) is the Green’s function of the first kind for the region 
under consideration. It is defined by the following conditions: 


~ 
~~ 
® 

i 
— 


(10.6) 


0°G 0°G 

— + — = 0 except at (&, n), 

Ox? dy? 

' (10.7) 
G(x, y; & 0) ~ — log {(x — &)? + (y — n)*}1/? near (€, n), 


G(x, y; &, ») = 0 over the solid boundaries. 
As is well-known, it has the reciprocity property 
G(x, y; & 0) = G(E, 0; x, 9). (10.8) 
For the case of a channel, it is given by the real part of 
F Tv Tv 
f~a=- {log sh a (2 — 20) — log ch z (z — iat, 


(g = «+ ty, 20 = E+ in, 2 = E — i). (10.9) 
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Let us now consider the behavior of a particular element of fluid at (£, 7) having 
an excess (or a defect) of vorticity corresponding to the secondary flow (10.5) and 
(10.6). It causes a distortion of the main vorticity distribution as indicated in Fig. 8. 
After a very small interval of time 6¢ the vorticity at the point (x, y) is changed by 
the amount 


d(x, y) = — v(x, y)btd (y), (10.10) 


because it is replaced by a fluid element from below, which retains its original vortic- 
ity. This change produces an effect at the “vortex,” i.e., at the element of fluid under 


(x,y) 


-_ 
ail 


oO x 
Fic. 8. Acceleration of vorticies in a non-uniform field of vorticity (¢¢ (y) >0, T'>0). 





consideration at (£, 7). It can be easily seen that the effect is a small velocity with 


components 


1 7) 

bu(t, n) = —{ f Gu, yi &, n)0o(x, y)dxdy, 
2a On 

(10.11) 


1 a 
bo(é, n) = — —{f ago ¥; &, )O¢(x, y)dxdy, 


the integrals being extended over the whole region between the planes. Dividing these 
quantities by 6¢ and passing to the limit 6-0, we have the following components of 


acceleration at the point (&, 7): 


1 rs) 
oe ee f f — G(s, 9: EWI, MEE Gedy, 
setsthaatrs (10. 12) 


1 fs) 
afi) = = f f 5g OCH 95 & MHC, INE (oDdady. 


Let us first consider the y-component of this acceleration. From the special form 
in which x and £ enter into the Green’s function [cf. (10.9)], we can also write 
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1 re) 
ay(é, 1) = —— f f vo'(x, ytd (y) —G(x, 9; &, n)dxdy. (10. 13) 
2 Ox 


If we multiply this equation by ¢’(, 7) and integrate over the whole region, we have 


the final formula 
ff xe mere asin = ff foe, 9} 206 (raady, (10. 14) 


upon using (10.6). Before discussing its signifiance, let us first notice that 


ff ote n)i'(&, n)dédn = 0 (10. 15) 


if ¢’(£, n) is an even function of §—&p, i.e., if the vorticity distribution ¢’(£, 7) has a 
symmetry about the line £=£. For then v’(x, y) is an odd function of x—£ and 
a-(£, 7) is an odd function of —£, both being the consequence of the fact that 
G(x, y; &, n) is an even function of x —£. Hence, we have the conclusion. 

If we recall that the vorticity ¢’(&, 7) is spread over a small region, we may take 
T=//f¢’(é, n)dédn as the strength of the “superposed vortex.” If we divide the left- 
hand side of (9.14) and (9.15) by I, we may consider the results as giving the compo- 
nents of the “average acceleration.” The x-component of acceleration vanishes; the 
sign of the y-component depends upon the sign of the superposed vortex and the sign 
of ¢/ (y). This component of acceleration is the one used in the above physical con- 
siderations. 

It should be mentioned that in considering the stability of a motion we deal with 
a vortex pair. Although this makes it difficult to obtain a compact formula for the 
average accelerations of the individual vortices, a kinematical consideration such as 
that given above (cf. Fig. 8) shows that the general tendency is not changed. Further- 
more, the two vortices are soon separated, because they are situated in layers of dif- 
ferent mean velocity. 

Another point should be mentioned. If we notice the tendency for the main vor- 
ticity to be swung around the secondary vortex, there is an acceleration of every ele- 
ment of fluid toward the vortex. Whatever this acceleration may be, it is expected 
to be of minor importance, because the effect is spread out over the whole field. This 
point will be brought out clearly in the following derivation of (10.14), where we shall 
study the whole phenomenon from the point of view of pressure forces. The accelera- 
tion will be identified with the negative of the pressure gradient divided by the density 
of the fluid, because the effect of viscosity has been neglected. Thus, if we can calcu- 
late the pressure disturbance corresponding to a given velocity disturbance, the left- 
hand side of (10.14) can be calculated. 

2) Second derivation, by consideration of pressure forces correlated with vorticity fluc- 
tuations. To calculate the pressure distribution from a given velocity distribution, we 
use for the pressure a differential equation of Poisson’s type obtained by taking the 
divergence of the equations of motion. Thus, if the equations of motion are’ 





10 The usual notation is used: x; (¢=1, 2, 3) are the coordinates, us are the components of velocity, 
p is the pressure, and p is the density of the fluid. Summation over a repeated index is understood. For a 


discussion of this type, see Lichtenstein’s book [26]. 
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Ou; Ou; 1 op 

















— + aw; = — — — + vAui;, it = 1, 2, 3), 10.16 
at * Ox; _ ( ( 
and the equation of continuity is 
du;/dx; = 0, (10.17) 
we have 
A(p/0) = - «, (10. 18) 
where 
Ou; OU; O( m1, u O(ue2, u O(u3, u 
ee in. — €:;€ij — Wifi; = { (11, U2) (u2, us) (us ae (10. 19) 
Ox; Ox; O( x1, X2) O( x2, x3) (x3, %1) 
and 
1 /du; du; 1 fou; Ou; 
tg <— + ), oj = — -— (10. 20) 
2 \0x; Ox; 2\0x; Ox; 


are the components of deformation and of rotation. If we integrate (10.18) under 
proper boundary conditions, o being known at the initial instant, we obtain the initial 
distribution of pressure. The initial acceleration field is then obtained from (10.16) 
as the negative gradient of the pressure, if we neglect the effect of viscosity. 

For a perfect fluid, the only boundary condition at a solid wall is 


umn; = 0, (10.21) 


where ; is the outward normal of the boundary surface. If we multiply (10.16) by 
n;, neglecting the effect of viscosity, we have 


1 op Ou; 
=V j (10.22) 


> tn ame o— &i, 
p On Os 
where Vo is the velocity along a stream line on the boundary, and ds is an element 


of its arc. If we write 
Ou; él; OVo 
ux = Voli, —= Ve— + |, —) 
Os Os Os 


where /; are the direction cosines of the velocity over the boundary surface, we have 


2 
1 ap V 
fae. (10. 23) 


where R is the radius of curvature of the stream line, R-!=n,01;/6s. This relation ex- 
presses the balance of pressure and centrifugal force. With a given distribution of 
velocity, the right-hand side is known. We have thus a potential problem of the second 


kind for the pressure. 
wo-dimensional flow between parallel solid walls. Returning to the problem at 


hand, we have the very simple boundary condition 


dp/dy=0 at y= +b. (10, 24) 
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Since the main motion is a two-dimensional parallel motion, we have 
= w(y) + u' (x, y), ue = v(x, y), u3 = 0, (10.25) 


where w(y) represents the main flow, and wu’ and v’ give a secondary flow approximat- 
ing a vortex. Equation (10.18) becomes 


42 2 , / , 
LY) Mn re eo 
We note that we can write o=01+02, where 
a flan ee .. (10.27) 
O(x, y) Ox 


a; depends upon the structure of the secondary vortex itself, and a2 depends upon its 
interaction with the main flow. We shall also separate the pressure into two parts and 
require them to satisfy (10.24) separately. Thus, 


p= pit ho, o = 01+ O2, 








1 /d? 0? fe) 

— (So -} =) = — 01, opr = 0 at y = + b, 

p \ dx? Oy? oy (10. 28) 
1 0? 2 0? 2 0 2 

= (+ *) = - 0, nae e 

p \ dx? dy? Oy 


We can reduce our problem to that of the first kind by looking for the acceleration 
a,(x, y) in the y-direction, ay= —(1/p) 0p/dy. If we differentiate (10.28) with respect 


to y, we have 








dy = a1 + a, ) 
07a 07a 001 07a 07a Oo. > (10.29) 
Ox? dy? 7 ay Ox? dy? fi oy J 


with a,=0, a2=0 at y= +b. The x-component of acceleration is zero, from symmetry 
considerations. The solutions of (10.29) are 


1 Oo}; 
oi(x, y) = — — | fe. y: &, o(S), aed 
1 O02 
x f foe. y; g, n) ($") dean, 


where the integrals are extended over the whole region between the planes. These 
formulae give the distribution of acceleration. Actually, it is more convenient to deal 


(10.30) 


ll 


a(x, y) 


with the integrated quantities 


I, = ff a(x, y)fo(x, y)dxdy, (10.31) 


I, = ff a2(x, y)fo(x, y)dady, (10.32) 
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JIi= ff a(x, y)f'(x, y)dxdy, (10.33) 


J2,= ff ao(x, y)f'(x, y)dxdy. (10.34) 


The first two integrals correspond to the accelerations of the main flow by the second- 
ary flow itself and by the interaction; the latter two quantities correspond to the ac- 
celerations of the secondary flow by itself and by the interaction. It can be verified 
(as will be done ape that 


--f fu- saad, (10.35) ot (10.36) 
a" (10.37) Se fe + asi (10.38) 


We note that (10.38) is essentially a reproduction of the formula sale the signifi- 
cance of which has been discussed above. The integral J; is equal to the negative 
of Jz. This is the above-mentioned acceleration distributed among the fluid elements 
throughout the field. It is therefore relatively unimportant. Thus, all the statements 
made in the last section have been verified, if we can verify (10.35)—(10.38). 

Verification of (10.35)—(10.38). To verify these equations, let us first examine the 
behavior of the quantities u’, v’, 0p/dx, 0p/dy for large values of x. From the expres- 
sion (10.9) for the Green’s function, we see that if ¢’ vanishes sufficiently rapidly as x 
becomes infinite, we have 


uw’ =O(|x|-*), vf = O(| x|-%), (10. 39) 
for large values of x. From the equations of motion, we then find that 
op op 
—=O(u') = O(| x|-), — =O’) =0(| «|-%. (10.40) 
Ox oy 


This will assure the convergence of the integrals involved and the validity of the steps 


taken in the following transformations. 
In the first method of derivation, we have been mainly concerned with Jz. We 
shall therefore consider it first. Referring to (10.30) and (10.5), we see that 


Oe2 
J» -ffve n) (>) didn. 
&” 


If we now introduce the value of o2 as given. by (10.27) and replace (&, 7) by (x, y), 


we have 
--2f five, = 


On integrating by parts with respect to x, we obtain 


J; = off v= (es o)dxdy -ff - (v'*&o) + 0”? Sh aay. 


The result (10.38) or (10.14) is thereby verified. 





“i fo)dxdy. 
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Similar calculations can be carried out for the integrals J;, J2, and J;. Thus, 


Oa, 
qT; = ff atutady = ff w(y) —dxdy 
oy 
1 0? 
-ff w(or+ — <2) andy, 
p Ox? 


by (10.28). If we note that 
d(u’, v’) ) Ou’ ) Ou’ dw dy’ 
w— = 2—([ wr’ —) — 2—[ wr’ —)- — — » 
O(x, y) Ox oy oy Ox dy dy? 


the above integral is easily transformed into the form (10.35). Following an exactly 


analogous process, we have 


1 O* De 
I,= ff w(o +—- ) axay = (), 
p Ox? 


when we make use of (10.27). The integral J; has also the significance that it is the 
effect of the solid boundaries upon a general flow y’(x, y) consistent with (10.40), 
because it is independent of w(y). Using (10.30), (10.5), and (10.27), we have 


00; 
Ji -ff aif’dxdy -{fve (=) didn 
OV t._ 


4 / , 


, (u,v) 
m=? w —————— daay. 
O( x, y) 
d(u’, v’) ) ( ) i) ( ~) 
; tv op 
s' — —a—i ee —ji~-—iser — Ih, 
O(x, y) oy Ox Ox oy 


we see that J;=0. The results (10.35)—(10.38) are thereby verified. We have thus 
completed the investigations indicated at the beginning of this section. 


If we note that 


NOTE ADDED IN PROOF. In a very early work, [Phil. Trans. Roy. Soc. London (A) 
215, 23-26 (1915)]| G. I. Taylor gave a physical interpretation of Rayleigh’s results 
on the stability of the laminar motion of an inviscid fluid, based on momentum con- 
siderations. He also indicated clearly that a motion, stable according to Rayleigh’s 
criterion, may be unstable through the effect of viscosity. 


(To be continued) 
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ON THE VIBRATIONS OF THE ROTATING RING* 


BY 


G. F. CARRIER 
Harvard University 


1. Introduction. An interesting addition to the group of problems dealing with 
thin elastic rings is the analysis of the vibration of a circular ring which is rotating 
with constant speed about its geometric 





axis. In this paper, the small bending vi- Ar T+2T de 
brations of the unconstrained ring are | nd 
analyzed and the frequencies at which [+56 ¢ 





such vibrations can occur are determined. 
For various problems of the partially con- 
strained ring, it is shown that the “free 
vibrations” differ essentially in character 
from those of the free ring, exhibiting a 
group of natural modes characterized by 
linear combinations of trigonometric func- 
tions. The forced vibrations of both the 
free and supported rings are also treated. 
2. The dynamic equations. The three 
equations needed to specify completely 
the plane motion of an element of a ring, 
such as the one shown in Fig. 1, are de- 
rived from a consideration of the forces 
and moments acting on the element and _Fyg. 1. Element of ring in initial and distorted 
the components of acceleration of the ele- positions. Center of rotation is 0. 
ment. The summation of forces along 
o’a’, the summation of moments about o’, and the summation of moments about a’, 
lead to this required set of equations, which is 


Or k k : 
— Rdd + |= + (p - <4) cos (a — v) — (v -—— Ar) sin (a — ») |ao =0, (1) 
c r Tr 


oR 1 0M k k 
- — +(¥ -— Ar) cos (a-9+(P — ~ 4.) sin (a — v) = 0, (2) 
00 Lal 00 r r 


ar 
R+ 56 de 


rr = — — rk — ——_ 0. (3) 


Here, the notations are as follows: @ is a polar coordinate of a point in the undeformed 
ring, referred to axes rotating with the ring; r and 7 are respectively the radii of 
curvature of the undistorted and distorted rings; 6, h and J are respectively the width, 
thickness and cross-sectional moment of inertia of the ring; R, r and M are respec- 





* Received March 2, 1945. 
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tively the tensile force, shearing force and bending moment, as shown; N/r and P/r 
are the components of the external forces in the directions of the tangent and normal 
in the undistorted state, as shown; A, and A» are the components of acceleration in 
the directions of the tangent and normal in the undistorted state, as shown; E is the 
elastic modulus; p is the density; k = pbhr?; c= Ebh; a= EI/r’. 

In Eqs. (1)—(3) the corrections, arising from the Poisson strains, for the moment 
of inertia expressions, etc., have been omitted as usual. The formulas needed to sup- 
plement the above equations are 


EI O°u 1 1 O7u 
M =— (we +u+- ), (4a) —=- (1 —-&m—-u-—--— =). (4b) 
r a6? ry r 00? 
Ov Ou 
R= Ebheg = clu tut ), (4c) a=—, (4d) 
06 00 
O7u Ov 
Ay = | 2 — wt + mw + 0) — 20], (4e) 
Oi? ot 
07y Ou 
Ae=r —+ 2w — — wp |, (4f) 
or? ot 
r Ov =0*u 
do = — (1 + e&)dd = (1 +—-- =) a, (4g) 
T1 06 06? 


‘where ruo is the radial displacement from the rest position to the rotating equilibrium 
position, 7u is the radial displacement from the rotating equilibrium position, rv is 
the tangential displacement relative to the rotating axes,' é is the tangential strain 
(€9=uo+u+0v/00), and w is the constant angular velocity of the ring. Equations (4a) 
and (4b) are the well known expressions for the bending moment and curvature, re- 
spectively, of a bent ring;? (4c) is a one-dimensional form of Hooke’s law; (4d) is the 
rotational displacement of the element, and is found by inspection of Fig. 1; (4e) and 
(4f) are the expressions for the radial and tangential components of acceleration, when 
u and v are referred to a rotating coordinate system;* and g is obtained from Fig. 1 
and Eq. (4b). 

The value of uo is obtained by writing u=v=0 in Eq. (1). After substitutions from 
Eqs. (4), it becomes Ebhuy = kw?(1+ 10), or 


uo = (1 + U0) kw?/c. (5) 


Throughout this analysis, we shall consider only those vibrations for which u and v 
are small compared to unity. We are therefore justified in disregarding terms in u’, 
uv, etc., as compared to u or v. In the limit, that is, as the amplitude of u and v tend 
to zero, the equations obtained in this manner would be exact. However, the equa- 
tions so obtained would still be encumbered by terms of the type uju, ugv, - --, in 
addition to those found below in Eqs. (6) and (7). We also neglect these terms since 








1 This rather unconventional notation is used to provide uo, « and v with dimensionless properties 
and thus produce somewhat less cumbersome equations. 

2S. P. Timoshenko, Strength of materials, D. Van Nostrand Co., New York, 1930, p. 459. 

* These are easily deduced from the vector forms given in L. Page, Introduction to theoretical physics, 
D. Van Nostrand Co., New York, 1941, p. 103. 
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they can produce no qualitative changes in the results and since they will drop out 
a yway when the procedure leading to Eq. (9a) is introduced. Finally, since the am- 
pli:udes of P and N must obviously vanish when u and v tend to zero, terms in 
Pu, Pv, - - - , must also be omitted. 

Following this procedure, using Eq. (3) to eliminate 7, and substituting from 
Eqs. (4) when necessary, we obtain from Eqs. (1) and (2) 


Ov 01 07 O74 Ov 
P-c¢ w+) — 0S + Se) + he —+u4- — 
06 064 06? 0g? 06? 


-E ov a 0% a < 








prc peers pat 
ot? ot c 368? c 36at? 


rs) ov O*u dou Ou 0*n Ou 
w+ 0< (ut =) —o( +) + bat = - #[ = + 20], (7) 
06 06 06" 00 00 ot? at 


We may easily arrive at a single equation in u only by performing on Eq. (6) the 


(6) 


operation L where, 

a? g2 . 
+k—») (8) 
06? Or? 





L=¢ 


solving Eq. (7) for L(v), and substituting the expression found by the latter step into 
that found by the former.‘ We utilize the abbreviations, 


s = th/a/k, p = wV/k/a, e = a/c = h?/12r?, 


and the equation resulting from the foregoing procedure takes the form, 


(2 -) +02 +2(2eieteog 
tL \ae2 as? 000s 36? \ae? 06? \ a6? 
a4 a2 / a (# 3? 
i. (5 * ‘ sae * 1) 2 “( i 3) as? 


on (+1 28) — + (ui » 2 —| 
+ —— ~ BP eee + OF = 9a = 
*\ 30? abas 36? 364 


as as ae 2a a a 
iS-nd-05-E-DeF 
: sac abas3 a6? «8* ~—sa02 as? 


(Ff 0/N oe (? 0/N 08 N 
LOC) EO mi Oba ° 

007 \ a 06 \a ds? \a Os \a 000s? \ a 
This equation governs the motion of the ring, provided the driving functions N and P 
do not imply that the’ deformations be large. The validity of this equation may be 
partially checked by considering a physically trivial problem. We consider the freely 
spinning ring (no supports) and suppose N to P to vanish. Under these conditions, 
the motion of the ring which is initially not deformed from its equilibrium shape is 
given by u=(a+8s) cos (@+us), v= —(a+8s) sin (8+ us). That is, the ring moves as 

















‘ A similar procedure will provide an analogous equation in v. 
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a rigid body with (dimensionless) angular velocity wu and translational velocity 8. 
Eq. (9) must and does allow this solution for all a and B. 

In each of those problems to be considered, Eq. (9), in its present form, leads to a 
solution which is simple im form but which requires the solution of unnecessarily long 
algebraic equations. This computational work may be eliminated at the expense of 
small errors in accuracy when we consider only rings for which h/r is small. In this 
case, terms in Eq. (9) of order € and e? may be neglected compared to those of order 
one. This leads to the following equation, which is exact for the limiting case where 
u, v, € tend to zero: 


_ a2 / 82 a? / 32 2 02 02 92 
Li(u) = | p?— ( — + 3) - — + 1) — 4u —— — | — — 1 }—_ lu 
06? \0@? 06? \a6@? 000s 06? Os? 
0? /P 0 /N 
--£0-20) “ 
06? \ a 06 \a 


The analogous equation in v is 


N 0 /P 
Li(v) = —+ =(—). (9b) 
a 06 \a 





In each of the problems to follow, the expressions obtained from Eqs. (9a) and (9b) 
for the natural frequencies and amplitudes are valid to within errors of order h?/r? 
for those frequencies of order (a/k)'/*. These vibrations may be termed “bending 
vibrations” since they are essentially inextensional forms of motion. When the fre- 
quencies are of order (c/k)'/, accurate results may be obtained by direct use of 
Eq. (9). In this paper, we shall use only Eq. (9a) or (9b) since all of the characteristics 
of the effects of rotation on the behavior of the ring will appear in the solutions so 
obtained. The one exception to this statement is found in connection with the dilatory 
vibrations, This type of motion can not be predicted by Eq. (9a) or (9b) because these 
equations are those for essentially inextensional motions.’ We shall, then, when in- 
vestigating this mode, refer to Eq. (9). For this mode, uw and v are independent of 0 
as must be those parts of P and N which excite such a motion. Hence Eq. (9) reduces 
for this case to 


04 c\ 0 0? /P 0 /N 
Lo(u) = | —-+ (su + — )=| “= ( -) + 2p <(—), (9c) 
os a / Os? Os* \ a Os \a 


or, in terms of 2, 
0? c N eis 
Lo(v) = (—+——-— ")— — 2y——-: (9d) 
ds? a a Os a 
3. The unconstrained ring. Investigation of the solutions of Eqs. (9a) or (9b) of 
the form 
u, = U, cos n(@ — B,s) + W, sin n(@ — 8,5) (10) 
yields, when N and P vanish identically, 


5 Equation (9a) may also be derived from Eqs. (6) and (7) by use of the assumption that u+0v/d0Ku; 
this excludes the dilatory motion. 
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Qu (n? — 1) eae oa 
‘. = + eerEE Eee aa /y2 2 1 =e n = ns 10a 
B aa rrr Ver + nv? + q p (10a) 


where p, and gq, are self-defining. 
If we further apply a set of initial conditions, such as 


Uy, 
u, = U, cos n0+ W, sin nO and _ = @Q, at s 
s 


0, 


to the two solutions defined for a given n by Eggs. (10) and (10a), namely 


tn = Un[an cos n(8 + gus — PS) + by cos n(0 + gus + prs) ] 
+ W,[c, sin n(0 + gus — pas) + dy sin n(O + gns + pas) |, (11) 


we obtain 


un = U, | cos n(@ + gns) cos np,s + = sin n(@ + q,S) sin nbes | 
p 


n 


x Yn . 
+ W, sin n(0 + nS) CoS nP,s — — cos n(6 + gnS) sin nas |. (12) 
Pn 
Each term in the foregoing bracket defines a possible free vibration of the ring which 
is unconstrained at all points against either radial or tangential displacement. Each 
of these terms may be interpreted as defining a “normal mode” of vibration, wherein 
a sinusoidal deformation of angular frequency mp, travels with respect to the rotating 
axes at an angular velocity’ —q,. The “nodal points” thus move with respect to co- 
ordinates fixed in the ring. The terms “normal mode” and “nodal point” have been 
used somewhat loosely here, but they adhere to the usual definitions of the terms if 
the motion is described relative to axes rotating with velocity 2, =w(m?—1)/(m?+1). 
For the stationary ring, the value of the angular frequency reduces to 


(a/k)*!*np, = n(n? — 1)[a/k(n? + 1)]/2, 


which is in agreement with previously derived results.” 

A solution to Eq. (9a) may be obtained for arbitrary initial distributions of radial 
deflection U» and radial velocity U@, provided these initial conditions do not imply 
an extensional motion. The restriction 


Qe Qe 
f U,(0)d0 = f Ug (0)d0 = 0 (13) 
0 0 


is certainly sufficient to insure this provision since it, together with the continuity 


requirements, 
2" Ov - oP 
f — do = f — =0 
0 06 0 000s 


allows u+0v/06 to vanish for all 6, at and shortly after time ¢=0. 
The restriction defined by Eq. (13) together with the requirement that u be con- 
tinuous, implies that Up and Ug may be expanded in Fourier series in which the 








® », and g, are, of course, dimensionless quantities which define the angular velocities. 
7 J. P. Den Hartog, Mechanical vibrations, McGraw-Hill, New York, 1940, p. 123. 
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constant term vanishes and in which no terms corresponding to »=1 will appear, 
since these terms define no distortion.* These series may be used in conjunction with 
series of solutions of the type given by Eq. (11). Thus the coefficients, and hence the 
motion, are determined. The motion due to rigid body displacements may, of course, 
be superimposed on such solutions. 

4. Forced vibrations. As will be seen in the following section, the investigation of 
the possible motions of the constrained ring requires, as a preliminary step, the 
determination of the behavior of a ring acted upon by a force distribution 
N/a=2A cos n@ cos ds. The problem involving a driving function P of similar 
form is obviously covered by this problem. If we split N into two parts, 


N/a = A[cos (nO — ds) + cos (n8 + ds) ], 
a particular solution of the form, 
v = b, cos (nO — Xs) + Cy, cos (nO + As) (14) 


is easily shown to exist by substituting this expression into Eq. (9b). The coefficients 
b, and c, are readily found when this is done, and are given by 


X,(A) = 6,/A = [n(n + + 1)(pra — q.) mo hon +1)” - 4nd], (14a) 
V,(A) = co/A = [n(n + 1)(p. — gn) — (2 +1) + 4rd], (14b) 
unless \ is one of the values given by \?=7263. B, is either of the values given by 


Eq. (10a). 
When N/a=2B cos n@ sin Xs, we have 


v = d, sin (n@ — ds) — e, sin (nO + ds), 
and 
d,,/B aed Ae e,/B —_ Y,, 

The quantities X, and Y, are useful later in the paper; hence the special notation. 

We see now that the motion of the unconstrained ring resulting from the type of 
loading described is composed of two waves of different amplitudes traveling around 
the ring with equal but opposite velocities. We note that there is, for each n, one value 
of \ for which there are fixed nodal cages in so far as tangential motion is concerned. 


This value of X is defined by X,+ Y, 
It follows from the linearity of our cavities that the driving function N of the 


more general form, 


N/a = ps Amn COS nO Sin Ams (15) 


will correspond to a solution 
= >> bmn Cos (08 — AmS) + Gmn COS (20 + AmS). (15a) 
Terms X,,, and Y,,, are defined as were X,(Am) and Y,(Am) in Eqs. (14a) and (14b). 
The particular problem in which the exciting force is given by 


N/a = Ao cos Xs (16) 





* This is seen in the discussion following the introduction of Eq. (9). 
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has no inextensional solutions. As was mentioned previously, we must, in this case, 

use Eq. (9d) for the determination of v. The solution has the form v=)» cos As where 

bo is given by 

bo X~? + €(1 — p?/d?) 

Xo(d) oo . 
Ago 1+ 3eu? — ed? 


(16a) 





The solution arising from the loading N/a=Ao sin As has the same coefficient. For 
small A?, v= A pA~? cos As. 

The natural frequency for the dilatory type of vibration is found by letting the 
denominator of Eq. (16a) vanish. Its value is given by 


hd = [(1 + 3eu?)/e]*/2. 


Returning momentarily to the question of accuracy, we note that here as in all 
subsequent problems the exact values of X, and Y, differ from those obtained in this 
section by terms (in the denominacor) of order ¢. Our work is accurate then when 
A\Keo/?, 

5. The supported ring. The first fact to observe in the investigation of the “free 
vibrations” of the partially constrained ring is that when WN and P vanish identically, 
no solutions to Eq. (9a) which obey the boundary conditions can exist. Specifically, 
we consider the ring to be supported by a number of evenly spaced, rigid, radial sup- 
ports (let there be J of them), and suppose the ring to be so fastened to these supports 
that radial motion is unconstrained at all points, but that v(27i/J, s) must vanish 
for all values of s and for each integer i. The first part of the appendix is devoted to 
the outline of a proof that Eq. (9b) has no solution under the foregoing conditions. 
Since the same proof holds for Eq. (9), we must conclude that the supports exert 
reactions which are to be accounted for in the differential equation by a function N 
which does not vanish identically. The problem, physically, becomes that of determin- 
ing what periodic forces, applied at the supports, are capable of sustaining a motion 
wherein the supported points of the ring have no displacement (tangentially) at any 
time. (We assume that the supports must move at precisely the speed w.) Mathe- 
matically, we must determine the eigenvalues \; and the corresponding solutions of 
the differential equation wherein we set P =0 and 


oo 


N/a = E +2 > cos nt] [A cos As + B sin ds]. (17) 
n=xJ ,2J,°+° 
This expression defines a loading which must correspond to a motion which has period 
2r/J in 0, since the force is the same at each support. When J is even, there may also 
be solutions periodic in +/J which don’t imply extensional motion. We shall not con- 
sider these, however, since both the procedures and results are analogous in the two 
cases. 
We have already shown [Eqgs. (15) and (15a)] that for loadings of the type given 
by Eq. (17), solutions of the form 
v(0, s) = >. [b, cos (n@ — As) +--+ — e, sin (nO — As) } (18) 


n 


exist for all \ except those for which \? =?8%. Using Eq. (9b) we determined all co- 
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efficients except those of index zero, and these were found with the aid of Eq. (9d). 
We may then, in this problem, replace Eq. (9) by the following: 

Lo(v’) = Lé (N/a), Li(v = v’) = (NV = N’)/a, (19) 
where the former equation is merely Eq. (9d), v’ and N’ are those parts of v and N 
which are independent of 6, and the latter equation is Eq. (9b). The solutions are now 


defined by Eq. (18) and the coefficients by Eqs. (14a), (14b), and (16a). 
If we now plot for a continuous range of X, 


o(d) = >) Xn(A) + Va(A) = DS (bn + en)/A = DS (dn + en)/B, 
n=J n n 
we find that the resulting graph (Fig. 2) contains two singularities corresponding to 
each n. Furthermore, there are two values of \ which we may associate with each n 
for which o(A) vanishes. One of these lies between the two singularities belonging 
to n; the other lies to the right of these values. We denote the smaller by X,, the 


larger by X,*. 





30 T 


| 
100 OCA) ] NX 
































3.0 
0 2.09 6.9 d 16.6 19.4 
Fic. 2. Response curve for ring driven at two points. o(A) is the amplitude of the motion of the 
points of application of the force. The broken section of the curve is drawn to the scale 1:5. For this 


curve, u.=3. 


We observe now that the motion at the points 27i/J is given by v(0, s) 
=:Ag(X) cos As+Ba(d) sin As. The values \; and \,* therefore define the frequencies 
for which the tangential displacements at the points of support are identically zero. 
They are therefore the desired eigenvalues. The motion of the ring for any \; or \* 
is given by 

0:(0, s) = >, AiXin cos (nO — dis) +--+ + BiYin sin (nO — Xs), (19a) 


where A; and B; are determined by the initial conditions. 

The question now arises as to whether linear combinations of the v; will always 
describe the motion arising from an initial set of conditions which are arbitrary except 
for the previously prescribed periodicity in 6. An outline of a proof that the v; are 
complete in this sense is included in the appendix. Since the natural modes of the 
possible vibrations are described accurately by Eq. (19a) only for the smaller values 
of \,, it may seem at first that this question of completeness is superfluous. However, 
the completeness of such sets of solutions provides an assurance that no other possi- 
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ble solutions of the equations have been overlooked in the analysis. In view of the 
fact that the foregoing procedure started with a guess as to the probable form of the 
support reactions, this indication that certain initial conditions would not lead to 
different types of motions is helpful. We conclude, then, that all vibrations, arising 
from initial conditions whose Fourier series expansions are such that the low index 
terms predominate, can be closely approximated by sums of the form, 


v(0, s) = p 5 Kw,(6, s). (20) 


We also conclude that the supports will continue to exert exactly those reactions re- 
quired to sustain this motion, or in other words, those forces required to prevent all 
tangential motion of the supported points of the ring. 

The problem of the radially constrained ring may be treated in a manner similar 
to the foregoing, with analogous results. When the ring is constrained at its supported 
points against both radial and tangential displacement, Eqs. (9a) and (9b) must both 
be used. Support reactions of the form 


N/a=A E + 2>° cos no | cosks, P/a= Bi + 2>> cos nt] sin As, 


and solutions of the form 
v= > a, cos (nO — As) + ag cos (nO + As), 
u = > b, sin (nO + rs) — b,! sin (nO — ds), 


are assumed to exist as before. This time we find four functions analogous to a(A) 
which enter the equations for the motion of the supported points. When this motion 
vanishes, these equations become Ao;(A)+Ba2(A) =0, Aas(A)+Box(A) =0. The criti- 
cal frequencies are defined by o1(A)o4(A) —a3(A)o2(A) =0. Since nothing essentially dif- 
ferent from the preceding results would be shown, the formulas for the terms in the o;, 
the explicit expressions for the v;, etc. are omitted. 

The forced vibration problem of the supported ring can now be easily treated. 
For example, let us consider the ring to be supported as in the first problem of this 
section of the paper, but to be loaded by a force distribution which may be expanded 
into the form N/a =)°.D, cos (n@+vs). The particular solution to Eq. (9b) corre- 
sponding to the loading N, is found as before, and the function 


v(0, s) = >> g, cos vs = G cos vs 


representing the displacement at the supports, has an easily evaluated amplitude G. 
Remembering that a support-reaction, 


N2/a=A E + 2>5 cos no] COS VS 


produces a motion at 0.=2zi/J which is given by 2(0, s)=Aa(v) cos vs, we may de- 
termine from the response curve (Fig. 2) the value of A such that Ao(v) = —G. The 
motion is then given by the solution to Eqs. (19) corresponding to the loading Ni+ N2. 
When more than one value of v enters the problem, the solution is changed only by 
the fact that the summation now takes place over two indices. 
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We note that when » in the problem just discussed is equal to one of the Xj, a 
resonant condition exists, as one would expect from the results of the preceding prob- 
lem. Also, we note that when the loading consists of a single term, cos m@ cos vs, and 
when + is one of the roots of X,,+ Y,,=0, the solution requires no support reactions. 

Perhaps the most interesting result of this analysis is the observation that, unlike 
most problems of this sort, the forces exerted by the rigid supports of the vibrating 
system must be included in the differential equation before the solution can be ob- 
tained. 

6. The elastically supported ring. A rather interesting eigenvalue problem arises 
when we consider the ring with elastic rather than rigid supports. Let us suppose 
again that the ring is unconstrained radially but that the supports resist the displace- 
ment of the points of attachment by a force, N/a = — Kv(0, s). 

Using Eqs. (19) as before, we find that the differential equations governing the 


motion now have the form, 


Lo(v’) = — KLé [v(0, s)], (21a) 
Li(v — v') = — K[v(0, s) [230 cos n6 |. (2ib) 
For solutions of the form, 
,= pe [an cos (v8 — Xs) +--+ — 6, sin (nO? — As) |, 


Eqs. (21a) and (21b) become 
a=- KFX,(X), ~=— KHX,,(A), 
Bn = KFY,,(A), bn = KH Y,(A), 
In these equations, F=)>.(an+B,), H=>>.(¥.+5,), and the X, and Y, are again 
given by Eqs. (14a), (14b) and (16a). 
When Eqs. (21c) are added by pairs and then summed over 2, the following re- 
sults are obtained: F[1+Ko(A)]=0, H[1+Ko(A)]=0. But if F and H vanish the 


solution is the trivial one; hence, 


(n = 0, J, 2J,---). (21c) 


o(k) = — 1/K. (22) 


This equation defines the eigenvalues and hence the natural modes at which the sys- 
tem may vibrate. We note that as K tends to infinity \ approaches that value found 
in the problem of the more strongly constrained ring (as it obviously should). As K 
tends to zero, the solution approaches that for the unconstrained ring. It is again easily 
shown that the set of eigen-functions obtained in this problem is complete in the pre- 
viously used sense. 

A final problem in forced vibrations follows easily from the foregoing. Let the ring 
be supported as above, but with the supports rotating at a speed w+ya sin As, 
where w is again constant. Briefly, we replace v(0, s) by v(0, s)+ cos As, on the right 
sides of Eqs. (21a) and (21b). The previously used procedures lead to the familiar 
set of solutions with the resonant frequencies obviously defined by Eq. (22). 


APPENDIX 


In this section, we wish to show, first, that Eq. (9b) has no solutions which are periodic in both 6 
and s and which vanish at @=kx/J, (k=0, 1,2,---, P=N20). 
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We assume a solution periodic in @ and write it as the sum of an even and an odd function (in @); 
v(0, s) = x(0, s) + (8, 5), 


where x(@, s)=x(—@, s) and y(@, s)= —y(—8, s). In order that v satisfy the specified requirements, both 
x and y must vanish at the prescribed points. 

In an abbreviated form we write Eq. (9b) as L2(v) =d*v/900s, where the operator L2 transforms even 
functions into even and odd into odd. This can be seen by inspection of Eq. (9b). The equation is now 
separable into two parts: 

0*y 0°x 


L(x) = (a) L(y) = 


’ 2 b 
00d obat (b) 








We operate on Eq. (b) with the operator ?( - --- )/80dt and substitute Eq. (a) into the result, obtaining 
o'x 


30701? (c) 


L(x) = 
which has solutions of the form, x =) n@n cos 0 exp tyns. All even, continuous, periodic, solutions of Eq. 
(c) may be written in this form but all of such solutions will fail to vanish at the specified points unless the 
a» vanish identically, since ¥m/yYn is irrational and cos 8 never vanishes at 6=0. Equation (a) then reduces 
to d?y/d0dt=0. It now becomes obvious that the solutions sought do not exist. 

Before showing that the functions derived as natural modes of vibration in the section on the sup- 
ported ring are capable of describing all motions periodic in 27/J which arise from arbitrary initial con- 
ditions, we introduce the following notation: 


Xin + Vin = din, s[— Xin + Vin] = Yin, Xin + Vin = Ging * (d) 


For any motion described by Eq. (20), the possible initial displacements and velocities may be written 


y :¥ [(A iPin + Aidin) cos "6 + (Bidin + Bi ¢in) sin nb |, 


v(@,0) = 
t=—0 neO,J,++* 
, (e) 
ov " 
— (0,0) = 2 2 [(Aidin + AT Vin) sin nO + (Badin + Bi yin) cos n6], 
os ‘ n 
where the A,, - - - , B;* are to be determined by the initial conditions. However, any initial conditions, 


periodic in 27/J can be written 


ov 
0(0, 0) = >- an cos nO + Bp sin nd, ro (6, 0) = Yn Cos nO + 5, sin n6. (f) 
s 


n 


This leads to the relations, 


mz (A idin + Al din) = Ging p 3 (Bidin + BY, in) = Bn, [se 


This set of equations may be considered as a group to be solved for the A;, - - - , B;* whenever such 
solutions exist. Except for special cases, such a set of equations always leads to a unique set of solutions 
corresponding to each set of an, +++, 5n. Since Eqs. (f) can express all the specified sets of initial condi- 
tions, we see that Eqs. (e) can also accomplish this purpose and hence the »; are complete in fthe sense de- 
fined above. 

The case where i equals zero requires a few additional words. The fact that there is no Oo or 59 seems 
to imply that we have two too few equations for the determination of the A,, - - - , B;*. Hewever, there 
is only one root of #(A) corresponding to m equal to zero. Henc e, the correct correspondenc between the 
ai, ++, 6:and the A;,--- , B,;* exists. 

The proofs outlined in this section are not claimed to be rigorous. They are presented merely to out- 
line the reasoning by which the two hypotheses might be proven if so desired. 
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NON-LINEAR THEORY OF CURVED ELASTIC SHEETS* 


E. BROMBERG anp J. J. STOKER 
New York University 


1. Theories of plane and curved shells which neglect bending.’ The problem to 
be treated here is that of determining, under certain conditions to be stated later, 
the stresses and strains in a thin curved elastic sheet in the form of a surface of revolu- 
tion held fixed at the edges and subjected to a uniform pressure normal to the sheet. 
The problem thus falls under the general class of problems treated in the theory of 
elastic shells. 

In order to obtain a theory of elastic shells which is manageable from the mathe- 
matical point of view, it is customary to make assumptions? of various kinds, in addi- 
tion to those of the general theory of elasticity. These additional assumptions are 
usually based on the hypothesis that the shell is very thin. In what follows we shall be 
interested in theories which result when the following assumptions are made: 

1) The strains due to the normal stress on elements parallel to the surface of 

the shell are small enough to be neglected safely. 

2) All stresses are constant over the thickness of the shell. 

The first assumption is almost always made by writers on the subject of thin shells. 
The second assumption of course rules out what are usually called bending stresses. 

A linear theory of shells, with a considerable number of practical applications, has 
been worked out on the basis of the above two assumptions. It is usually referred to 
as the membrane theory of shells.* The salient feature of the theory is that it is “stati- 
cally determinate” since the stresses can be obtained from the equilibrium conditions 
alone without reference to the elastic deformations. This results in a very great sim- 
plification, by comparison with theories which do not neglect bending. However, the 
simplification is coupled with at least one rather serious disadvantage: it turns out 
that it is not possible to satisfy the kind of boundary conditions which it would be 
natural to impose in these problems, since the order of the system of differential equa- 
tions is too low. For example, the condition of a fixed edge (that is, the condition 
requiring the displacements at the boundary to vanish) cannot be satisfied in general. 

Most writers on the membrane theory of shells attribute the difficulty regarding 


* Received May 1, 1945. 

1 The theory developed in this paper is an outgrowth of a research project carried out by the College 
of Engineering of New York University, under a contract with the War Production Board. The investiga- 
tion, which was largely experimental in character, was concerned with the feasibility of constructing 
buildings circular in form with a thin steel roof supported by excess air pressure on the inside of the build- 
ing. The design problems which arose led to the theory presented here. In this case the sheets considered 
were so thin that there was no doubt about the validity of neglecting bending stresses. 

2 For an exhaustive classification of the very numerous possibilities here, see the recent paper of 
Chien [3]. (Here and in what follows, numbers in square brackets refer to the bibliography at the end of 
the paper.) 

* For full treatments of this theory and references to the literature, see the books of Fliigge [5] and 
Timoshenko [12]. 
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boundary conditions to the fact that bending is neglected, and it is true that use of 
the linear bending theory does make it possible to impose physically reasonable 
boundary conditions. However, there are cases in which the shells are so thin that 
the bending stresses are small compared with the “membrane” stresses.‘ It seems not 
to have been noticed that a theory which neglects bending stresses, but which neverthe- 
less makes it possible to satisfy physically reasonable boundary conditions, can be obtained 
by taking account of certain non-linear terms in the relations for the strains as functions 
of the displacements. This paper has as its main purpose the development of such a 
non-linear theory. 

Our theory is a generalization of an already existing non-linear theory for the case 
of a plane sheet® supported in some way at its boundary and subjected to normal pres- 
sure p. It is useful for our purposes to discuss the theory of plane sheets from a num- 

% 


Tro 


fe 
pet 
Txy 





Ss 


Tey 



































Fic. 1. 


ber of different points of view, with the object of comparing and contrasting this 
theory with the theory of curved sheets to be presented later. The undeformed posi- 
tion of the sheet is taken as the xy-plane, the system of stresses in the sheet is denoted 
by o:, gy, and tz, and the displacement components by u, v, and w. In Fig. 1 the 
notation for the stresses ¢, and og in polar coordinates (7, @) is also indicated. 
The equilibrium conditions for the stresses oz, gy, and Tz, in the sheet are 
O¢s OFsy OTzy  O6y 


—- = Q, - +— = 0. (1.1) 
Ox oy Ox oy 











The equation of equilibrium for the direction normal to the sheet is 


0°w 0°w 0*w 
Cy —— t-2r sy ho = p/h, (1.2) 
Ox? Oxdy dy? 
‘ The present investigation was prompted by the necessity of dealing with just such a case. 
5 The word “sheet” is employed here in a noncommittal way. In the course of our discussion a more 
precise significance will be given to the phrase “theory of thin sheets.” 
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where h is the thickness of the sheet and p is the pressure. The non-linear character 
of the theory under discussion stems from the retention of certain quadratic terms in 
the relations for the strains €,, €y, Yzy in terms of the displacements * 


Ou 1 /ow \? Ov 1 /dw \? 
oy == + ne sf €y =e a — = as 
Ox 2\0x oy 2 \dey 
Ou Ov Ow Ow 
Ts * + ee 


Oy Ox Ox dy 


(1.3) 


The stresses and strains are assumed, as usual, to obey the stress-strain relations 


1 1 2(1 + ») 
és = — (oz — Voy), a (cy — voz), ae“ a aes 


E E E 





where E and v are the modulus of elasticity and the Poisson ratio, respectively. 
Finally, we have the “compatibility” equation 


” f 0’w\? dw d*w 
Vo.+to, =E {( ) — = , (1.4) 


Oxdy Ox? dy? 





in which V? is the Laplace operator.? This equation is an integrability condition for 
Eqs. (1.3), expressed in terms of o:, ¢,, and w by the use of (1.1) and the stress-strain 
relations. 

Two different kinds of conditions will be considered at the boundary C of the 
sheet. In one case we prescribe the displacements #@, 3, w at the boundary; 


ui = u(C), d= r(C), w= 0. (1.5) 


In the other case, instead of the displacements # and i we prescribe the normal and 
shear stresses ¢, and 7, at the boundary; 


fn = g,(C), 7, = 7,(C), w= 0. (1.6) 


That we may impose the boundary conditions (1.6) is clear; the differential equations 
(1.1), (1.2), and (1.4) together with the boundary conditions (1.6) constitute the 
complete formulation of a boundary value problem for the determination of the 
functions o2, ¢,, Tz, and w. That the conditions (1.5) may be imposed could be seen 
readily by formulating our problem in terms of the displacements u, v, w alone, but 
we refrain from doing so here. The conditions (1.5) mean that the edge of the sheet 
is stretched in its plane by a fixed amount, which does not depend upon the applied 
normal pressure p. The conditions (1.6), on the other hand, mean that the stress at 
the edge is held fixed while the displacements there will depend upon p. 





6 Notice that of the three quadratic terms occuring in the usual expression for the strains, only the 
one involving w is retained in (1.3). The motivation for this is that the order of magnitude of the dis- 
placement w normal to the sheet can be expected to differ from that of the displacement parallel to the 
plane of the sheet. The experimental results (see the paper by Eck [4]) confirm the validity of this assump- 
tion from the physical point of view. 

7 These differential equations were first obtained by Fépp! [6] in 1907. They can also be obtained by 
neglecting the terms referring to bending in the non-linear theory of plates developed by v. Karman [11]. 
The equations have been solved by Hencky for the case of a circular sheet [9] and a rectangular sheet 
[10]. Bourgin [2] has treated the case of the rectangular sheet by methods different from those of Hencky, 
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For the purpose of comparison with our new theory of curved sheets, which we 
develop only for the rotationally symmetric case, it is convenient for us to formulate 
the plane sheet theory in polar coordinates assuming all quantities to depend only on 
the coordinate r. 

The differential equations corresponding to (1.1), (1.2), and (1.4) are for this case 


d 
=~ (ro,) = 4, (1.1)’ 
1 d (ro =) ar p/h (1.2) 
r dr dr : ; 
d E {dw\? 
sie a (ro4) +o, = (=). (1.4)’ 


For the purpose we have in view it is not necessary to write down the strain-displace- 
ment and stress-strain relations in polar coordinates. The boundary conditions (1.5) 
at the edge r= R become 

a = u(R), w = 0, (1.5)’ 
in which @ refers to the radial displacement at the edge. The alternate boundary con- 
ditions (1.6) become 


6, = ,(R), w= 0. (1.6)’ 


We consider three different specializations of the non-linear plane sheet theory as 
a basis for comparison with the theory of curved sheets to be developed later. These 
are: Case (a), a direct linearization of the differential equations; Case (0), the classical 
linear membrane theory; Case (c), the problem of Féppl-Hencky. We proceed to dis- 
cuss these three cases in order. 

Case (a). A direct linearization of the differential equations. If we simply neglect 
the non-linear terms in (1.2) and (1.3) we obtain the relations V?(¢.+¢,)=0, p=0. 
The sheet is therefore not deflected laterally; it is simply in a state of plane stress. 
From our point of view, such a linearization thus leads to a “trivial” problem. 

It is worth while to point out that the solutions for Case (a) are also solutions of 
the non-linear sheet theory if we impose the condition that the normal pressure p be 
everywhere zero. 

CAsE (b). The classical linear membrane theory. The well-known linear theory of 
tightly stretched plane membranes can be obtained from the non-linear sheet theory 
as an approximation to the solution of the boundary value problem in a special case. 
The approximation, as we shall see, results from a development in the neighborhood 
of Case (a). The special case of the non-linear theory in question arises when the 
boundary condition is taken in the form (1.6) with ¢, assumed to be a constant ¢>0, 
7, to be zero: * 

¢,=a>0, tT = 0, w = 0. (1.7) 
Furthermore, we make the important additional assumption that the applied normal 
pressure p is small compared with ¢. In other words we assume the membrane to be 
tightly stretched and then deflected by a relatively small normal pressure. 

We can solve this boundary value problem by a perturbation method consisting 
of a development in the neighborhood of the solution for the case in which w=0, 
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p=0 (that is, in the neighborhood of the undeflected state of the stretched sheet). 
The well-known linear membrane theory results as the second step in such a develop- 
ment. We need only develop az, oy, Tz, and w in terms of the parameter e defined by 


e=p/c, as follows: 


(0) 2 (2) (0) 2 (2) 
Cz=0, teéeoz t:°°, Oy=Cy téeoy +°°::, 

(0) | 2 (2) a), 8 (2) a (1.8) 
Tey = Try TETey Fe", w=ew tew +°°-, p = eG. 


The stresses (including the stress ¢ at the boundary) are of lower order in € than the 
deflection w and the applied pressure p. Insertion of relations (1.8) in the differential 
equations (1.1), (1.2), (1.4) and the boundary conditions (1.7) leads to a sequence of 
linear boundary value problems for the determination of the coefficients in the per- 
turbation series. For the terms of zero order in the stresses one finds readily the solu- 
tion o =o =¢, 7%) =0; in other words the zero order terms represent a state of 
uniform tension throughout the sheet. The zero order terms are also, evidently, the 
solution for the linearized sheet theory of Case (a). The differential equation for w 


is then readily found to be 
V2w) = — 1/h, (1.9) 


while the boundary condition is, of course, 
wh(C) = 0. (1.10) 


Equations (1.9) and (1.10) are those of the classical linear membrane theory (for 
unit normal pressure). For the applicability of this theory the essential condition is 
that the applied pressure p should be small compared with the initial stress ¢ in the 
sheet. We note also that this theory results when the stress is prescribed at the bound- 
ary rather than the displacement in the plane of the sheet; in other words, the linear 
membrane theory requires that the edge of the sheet be free to movée in the xy-plane. 

CasE (c). The problem of Fippl-Hencky. The boundary value problem which leads 
to our Case (c) is that resulting from the choice of (1.5) as boundary conditions for 
the non-linear sheet theory. This theory is sometimes referred to as the large deflection 
theory of membranes. It is not assumed, as in the above Case (b), that the normal 
pressure ~ is small compared with the initial stress in the sheet. In fact, we assume for 
the Case (c) that the displacements # and é at the boundary as well as w# are zero. 
We shall refer to this problem® on occasion as the problem of Féppl-Hencky. Our 
boundary conditions of course mean that the sheet was initially unstrained. Thus the 
stresses in the sheet are built up only as the normal pressure p is applied, and conse- 
quently the procedure outlined above for Case (b) is entirely inapplicable. 

As already stated, our purpose is to generalize the non-linear sheet theory (c) to the 
case of curved sheets. The essential step for this purpose consists in developing suitable 
non-linear strain relations for the curved sheet analogous to those (cf. (1.3)) for the 
plane sheet. However it is not entirely clear a priori in the case of curved sheets just 
which of the quadratic terms in the strain equations should be retained and which 
rejected. Section 2 is devoted to a derivation and discussion of the strain expressions 





8 The previously cited papers of Hencky [9, 10] and Bourgin [2] are concerned with this problem. 
The problems for the case in which the boundary displacements @ and @ are not zero (i.e. the case of an 
initially stretched sheet) appear not to have been treated. 
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used later as the basis for our theory. The discussion is confined to the case of the 
rotationally symmetric deformation of a surface of revolution. Only two displace- 
ments are involved in this case, the displacement u along a meridian and the displace- 
ment w along the normal to the sheet. 

Once expressions for the strains in terms of the displacements are available, 
it becomes possible to set up the integral for the potential energy in the sheet (as- 
suming Hooke’s law to hold) in terms of the displacements u and w. The equilibrium 
conditions can then be found as the Euler variational equations minimizing the po- 
tential energy. The result is a pair of second order non-linear differential equations for 
u and w which permit the boundary condition 4=#=0 for a fixed edge to be imposed. 

For most purposes it is, however, more convenient to formulate the curved sheet 
theory in terms of the two components o» and a, of the stress in the sheet along and 
perpendicular to a meridian curve, respectively, and the displacement w normal to 
the sheet, rather than in terms of the two displacements u and w. This is carried out 
in Section 3. In Section 4 the general theory is specialized for the case of the spherical 
sheet. The result is a set of differential equations for the curved sheet analogous to 
(1.1) and (1.2) for the plane sheet. In Section 5 a simplification in the theory for the 
spherical sheet is introduced which is valid for a spherical segment of small curvature 
(and probably also for all cases of spherical sheets). The differential equations of Sec- 
tion 5 are 

d 


— (ogsin@) = o, cos 8, 
dé 


1 d _ dw Rp : 
—( os sin 0 —]} = — {— + og + oe) sin 8, 
R dé dé h 
! (1.11) 


d 
— — (o¢ tan 6) + (1 + »v tan? Ooze 


a 
E { 294 , dw Pe 1 (<= )} 
= — 4w tan ae ep ee : 
R d6 2R\ dé 





The independent variable @ is the latitude angle measured from the pole of the sphere. 
These equations are exactly analogous to Eqs. (1.1)’, (1.2)’, and (1.4)’.® We refrain 
from writing the stress-strain and strain-displacement relations which are needed for 
a complete formulation of the problem. 

As boundary conditions at 6 =6) we assume either 


u(8o), w = w(4o), (1.12) 


u 
or 
Ge = 000), w= w(). (1.13) 
At the pole @=0 we require all quantities to remain finite. 

We wish to consider the three special cases in connection with Eqs. (1.11) which 
are analogous to the three cases discussed above in connection with the plane sheet. 
These are: Case (a), a direct linearization; Case (b), the analogue for curved sheets of 





® If we were to allow R to tend to infinity while @ tends to zero in (1.11) in such a way that Ré->r, 
the result would be the differential equations (1.1)’, (1.2)’ and (1.4)’. (The normal pressure p and the dis- 
placement w are taken positive in the direction toward the center of the sphere.) 
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the classical linear membrane theory; Case (c), the analogue of the Féppl-Hencky 
problem for curved sheets. We consider these cases in order. 

CAsE (a). A direct linearization. In contrast to the corresponding case in the plane 
sheet theory, we observe that neglect of the non-linear terms in (1.11) does not lead 
to a trivial problem. We obtain, in fact, the equations 





, 


d : Rp 
— (ce sin 0) = ao, cos 8, Coton, = — 

dé h 
(1.11)’ 


= (a, tan 6) + (1 +» tan®6) —( a+ o—) 
a , tan tan* = —I| w tan’ - tan 0— e 
ao , er d0 


These are the differential equations of what is called the “membrane theory” of thin 
shells. One observes that the order of the system (1.11)’ is two less than that of (1.11) 
The stresses can be obtained from the first two equations of (1.11)’ without reference 
to the strains and displacements. It is, however, not possible to satisfy in general the 
kinds of boundary conditions which would be natural in the physical situations en- 
countered in the applications. One such case is that of a fixed edge, which would imply 
the condition w=0 at the boundary. That this condition cannot be satisfied, at least 
in the case =const., is readily seen: the only solution of (1.11)’ that is finite at @=0 
is given by 

7=0%=>-—») w= pR*(1 — v)/2Eh. (1.14) 


All three quantities are constant throughout the shell. 

As in the corresponding Case (a) for the plane sheet, the solution (1.14) of the 
linear equations (1.11)’ is also a solution of the non-linear equations (1.11) if proper 
restrictions are imposed. Instead of prescribing the boundary values @¢ and @ in (1.13) 
arbitrarily, we would require them to have values consistent with (1.14). 

CAsE (b). The analogue for curved sheets of the classical linear membrane theory. The 
theory of curved sheets analogous to the classical linear membrane theory for plane 
sheets seems not to have been developed. For the case of a spherical sheet we can ob- 
tain such a theory from Eqs. (1.11) with the boundary condition (1.13), in which, 
however, we assume @¢ and W to have values consistent with (1.14) and set p=p. 
However, we assume that the pressure p in (1.11) is given by 


p= pO + ep, (1.15) 


in which ¢€ is a small (and, of course, dimensionless) parameter. The theory we desire 
then results from the terms of first order in the development of the solution by per- 
turbations with respect to e. We set 


(0) (1) 2 (2) 
w=w +ew +ew meee 


0 1) 
oo = 09 +05 a j (1.16) 
eh Pe 


and insert these series together with (1.15) in the differential equations (1.11) and 
the boundary conditions 


Go = 6g = — P™R/2h, w= pOR*(1 — v)/2Eh. (1.17) 
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The terms of zero order in (1.16) are then readily found to be those which would re- 
sult from the linear theory—in other words of”, oD, and w” have throughout the 


sheet the constant values prescribed at the boundary in (1.17). The first order terms 
are then found to satisfy the linear differential equations 


d 
= fos sin 0) = os. cos 0, 


1 d (0) dw) (“E (1) ’) 
-- sin 6 =— sin 0, 1.18 
“(os 70 ) h +oo + T¢ , ( ) 








d (1) 2 (1) E dw") (1) 2 
——(o, tané)+(1+rtan 0)oo = —{ tané +w tan 6}, 
10 R d@ 


¢ 








é 


and the boundary conditions 


6 =0, w =0. (1.19) 


Equations (1.18) and (1.19) are analogous to (1.9) and (1.10) for the corresponding 
case of the plane sheet. It must, however, be admitted that this “theory of tightly 
stretched membranes” for the sphere is somewhat artificial because of the fact that 
the “stretched” state is one for which the initial radial displacement w cannot be held 
zero at the boundary. 

CAsE (c). The analogue of the Fippl-Hencky problem for curved sheets. The differ- 
ential equations (1.11) are to be solved for a prescribed pressure p when the edge of 
the sheet is considered fixed, i.e., under the boundary conditions #=0 and #=0. 
In this particular case the condition #=0 can be replaced by the condition that the 
strain in the direction of the boundary curve is zero, which implies the condition 
4 —vee=0 on oy and a» at the boundary. The analogy with the corresponding case 
for the plane sheet is, as we see, exact in every respect. 

It should now be apparent that some such term as “sheet theory” is needed in 
addition to the term “membrane theory.” This is brought out by Table I which lists 
the Cases (a), (b), and (c) together with the present terminology. As we note, the 
phrase “membrane theory” is already applied to cases which have almost nothing in 
common. Consequently we would recommend (in accordance with a suggestion made 
by Bourgin [2]) that all of these theories which neglect bending be referred to in 

















TABLE I, 
tm Plane | Curved 
(a) Plane stress | Membrane theory of shells 
= (b) Membrane with small deflections | — 
a (c) Large deflection theory of membranes | — 





general as sheet theories. The Cases (a) and (c) could then be referred to as linear and 
non-linear sheet theories respectively, while the term membrane theory might be re- 
served for the Cases (b), i.e., for theories of initially stretched sheets which start with 
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the linear sheet theory as a first approximation and then proceed to a second approxi- 
mation by a development in the neighborhood of the solution to the linearized prob- 
lem. This terminology will be used in the remainder of this paper. 

In Section 3, the differential equations for the curved sheet theory are obtained, 
as we have already indicated, by variational methods. In that section also, the sta- 
bility of the extremal solutions for both the linear and the non-linear sheet theories 
[Cases (a) and (c)] is considered. At first sight one would be inclined to think that 
the solutions in the two cases would not differ greatly as far as stability is concerned 
if the pressure, thickness, etc., are the same in both cases. This is, however, not true. 
On the basis of the linear curved sheet theory, the solutions would appear to be stable 
whether the pressure / is inward or outward, that is, whether the sheet is in compres- 
sion or tension, respectively. On the basis of the non-linear curved sheet theory, how- 
ever, the solutions are unstable when the pressure is such as to cause the stress a 
in the sheet to be a compression.’® This result follows through consideration of the 
Legendre condition for our variational problem. In the case of the spherical sheet the 
stress g», as given by the linear theory, is a compression when the normal pressure p 
is positive (i.e., when the pressure is directed toward the center of the sphere). It also 
seems certain that the non-linear theory will yield the same relation between the signs 
of a, and p for the case of the spherical sheet, unless the displacements are very 
large. Consequently, we have assumed in our numerical work that the pressure p is 
negative, i.e., is directed outward, in order to avoid unstable cases. 

In Section 5 the non-linear curved sheet theory [Case (c) ] is formulated in detail 
for the special case of a spherical segment of small curvature. The differential equa- 
tions for the spherical sheet can be solved by power series in the independent variable. 
Graphs showing the distribution of the stresses and the normal deflection w along a 
meridian in a particular numerical case are given in Section 6. Perhaps the most strik- 
ing feature of these results is that the non-linear sheet theory [Case (c) | yields results 
which do not differ greatly from those of the linear theory [Case (a) ] except near the 
edge of the sheet. In particular, the stresses and the normal displacement w are nearly 
constant over most of the interior of the sheet, but change rather rapidly near its edge. 

This observation indicates that we have to deal here with a boundary layer effect. 
In Section (7) the existence of such an effect is deduced and treated explicitly. It 
turns out upon introduction of proper dimensionless variables in the original differ- 
ential equations that only one parameter x remains in the transformed differential 
equations. The quantity x is given by 


x = pR/Eh, (1.20) 


in which # is the normal pressure on the sheet, R the radius, E the modulus of elastic- 
ity, and h the thickness of the sheet. If x is allowed to tend to zero in the transformed 
differential equations the result is in the limit the differential equations of the linear 
sheet theory with a consequent lowering of the order of the system. Hence some 
boundary condition must be lost at the edge on the transition to the value x=0. The 
solutions of the boundary value problem for x0 can therefore not be expected to 
converge uniformly at the boundary to the solution of the problem for x=0. It is 

10 Comparison with the analogous cases (a) and (c) for the plane sheet theory is illuminating. It is 


clear that the plane sheet would be stable under edge compression if no lateral deflection were to be per- 
mitted, but decidedly unstable under compression if such a lateral constraint were not imposed. 
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possible to treat the boundary layer phenomenon by introducing a new independent 
variable which depends upon x in such a way as to stretch the boundary layer to 
infinity as x0, with the result that the convergence of the solutions becomes uniform 
with respect to «x at the edge. One notices that the value x =0 corresponds, according 
to (1.20), to the value zero for the pressure p. 

The boundary layer solution is given in this case very simply by an exponential 
function. It could be used to estimate the stresses in practice in cases for which x is 
small (and, of course, negative). In the usual cases it is not difficult to see that «x will 
be of the order of —0.0005 in practice, since pR/2h is the stress when p is constant, 
according to the linear theory, and hence x is a quantity of the order of the longitudi- 
nal strains." 

It is clear that the non-linear sheet theory could be worked out in detail rather 
readily in other cases such as those of the cylindrical and conical sheets. It would 
also be of interest to consider the case of the spherical 
sheet with a hole, so that two distinct boundary curves Y 
would exist. Various combinations of boundary condi- 
tions at the two edges should be considered; boundary 
layer effects could then occur at both edges. Cc 

From the point of view of the practical applications, 
another question is of interest. It is clear that bending 
effects will dominate the “sheet effects” near the edge of 
the sheet if the sheet is thick enough. This question is [ \//“w 
under investigation at the present time. / 

2. Expressions for the longitudinal strains. We as- 
sume the curved sheet to be the surface of revolution ob- 
tained by rotating about the y-axis the meridian curve C, Fic. 2. 


x= a(t), y= y(8). 6-9 


The parameter £ is taken to be the arc length of the curve. We consider only deforma- 
tions which preserve rotational symmetry, so that the deformation is completely de- 
scribed by the displacement components u and w along the meridian and along the 
normal to the surface respectively. 

It is convenient to introduce the angle 6 between the y-axis and the normal to the 
meridian. These notations are indicated in Fig. 2. 

The longitudinal strains in the sheet are defined in the usual way. If ds' is the de- 
formed length of the line element originally of length ds, then the strain ¢ in the direc- 
tion of the element ds is defined by 


(=) - 1+ 2¢. (2.2) 


ds 








It is useful to introduce the following relations between the original position (x, y, 2) 
and the deformed position (x', y', s') of any point P on the sheet: 

it The effect of the edge constraint seems to be such as to cause the stresses at the edge to be lower in 
value than in the interior of the sheet. Thus it seems likely that the usual practice in engineering design 
of ignoring the edge effect leads to estimates for the stresses which are too high, i.e., are on the side of 
safety. Of course, we are entitled to draw this conclusion here only in case the sheet is in tension. 
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x1'=x2+ucosé— wsin8, yi = y— usinéd + wos 8, zs} = ¢g. 


By making use of these relations the strains € and ¢€, in the direction of a meridian 
and a parallel (i.e., a curve 6=const.) are easily computed by using (2.2). We obtain 








1 /du 1 dw\? 1 1 /du 2 
é6 = —(|— — w)+ -| u-+—)+— —(— — w}, (2.3) 
p \dé 2p? d@ 2 p* \d0 
u cos 6 — wsin 6 1 /ucos@ — wsin6\2 
€g = —————_ + — : (2.4) 
x 2 x 


The quantity p is the radius of curvature of the meridian curve and x is, of course, 
the abscissa of the point P. 

Just as is done in the analogous case of the plane sheet, we retain only certain of 
the quadratic terms in the strain expressions, which then amounts to the assumption 
that these non-linear terms are considered to be of the same order as the linear terms. 
Thus it would be logical to reject the third term on the right hand side of (2.3) and 
the second term on the right hand side of (2.4), since they are squares of the linear 
terms. We shall follow this procedure and thus take for the strains the expressions 








1 (= )+ 1 ( A ), (2.5) 
é=—|— — w u+—}, : 
"pp \d0 2p? d0 

ucosé — wsin# 
€, = : (2.6) 

x 


The following special cases are of interest: 
a) The sphere. Here p=R (the radius of the sphere), and we find from (2.5) and 
(2.6) that 


es ee a — 

R \d@ 2R?\ dé 

é,= hs (u cot@— w). (2.8) 
R 


b) The circular cylinder. Here p= ~, 0=2/2, x=a (the radius of the cylinder), 
and &=y. We find that 


os (=). (2.9) es (2.10) 
&=— —|— }, ; &=--—-: 4 
‘ dy 2\Xdy . 


c) The circular cone. Here p= ~, 6=y=const. The strains are given by 


du 1 dw \? 
co = cos y+ — cost y(— : (2.11) 


dx 2 dx 


u cosy — wsin y 


és = 





x 





12 These expressions coincide with those used by Friedrichs [7]. Similar expressions were used earlier 


by Biezeno [1]. 
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3. Formulation of the boundary value problem for a sheet with fixed edges. On 
the assumption that Hooke’s law holds for the relation between the strains’ €s and €, 
defined above and the corresponding normal stresses o» and o4,"* we have for the po- 
tential energy V of the sheet the expression 


o&/1 Eh 

V= an f (> ——— [G+ & + 2veve,| — pu) xdé, (3.1) 
: \2 1-—¥» 

in which p is the normal pressure on the sheet counted positive in the direction of the 

inward normal. The quantities and x(£) are the arc length and abscissa of the me- 

ridian curve, and h is the thickness of the sheet. The quantities v and E are the elastic 

constants. 

The potential energy could be expressed in terms of the displacements u and w 
by replacing €» and ¢€, in terms of these quantities through (2.5) and (2.6). The varia- 
tional equations for the minimum problem would then clearly be a pair of non-linear 
ordinary differential equations for u and w, each of which would be of the second 
order. We shall not write these equations down since in the following we wish to work 
with the stresses a, and oy, and the displacement w as dependent variables. However, 
we do wish to draw one conclusion from the existence of two such equations. The 
differential equations for u and w are of the proper order to permit imposition of the 
boundary condition u=w=0 appropriate for a fixed edge. 

The variational equations resulting from (3.1) are 


, «a (x0) [“(+ eed (3.2) 
ms dé XO@ , dé , ao — ° 
; dw u un) o sin 0 p 
z o| 20o( = + *)]4 [24 "4 A] =o (3.3) 
x dé dé p p x h 


The quantity p in these equations represents the radius of curvature of the meridian 
curve; the quantities p, @, and x are, of course, given functions of £. In deriving (3.2) 
and (3 3) use was made of the stress-strain relations 


Ees = 06 — v0¢, Ee = 04 — voo, (3.4) 


and of (2.5) and (2.6) in order to introduce a» and o, as dependent variables. Equa- 
tions (3.2), (3.3), (3.4), (2.5), and (2.6), together with appropriate boundary condi- 
tions, yield the complete formulation of the boundary value problems we consider 
here. We note that there are six equations for the six quantities u, w, Oo, 0%, €0, €4- 

For the most part, we are concerned with the case of a sheet without a hole at the 
axis of symmetry, so that the quantity x in (3.2) and (3.3) has the value zero where 
the meridian curve crosses the axis, which we may assume to occur for §=0. In this 
case we would require the solution to be regular at £=0. At an edge =p of the sheet 
we require u=w=0, for a fixed edge. In view of (2.6) we see that this implies e,=0; 
hence we may prescribe the following conditions at a fixed edge: 


w= 0, (3.5) 


ss a, (3.6) 


In this way we express the boundary condition in terms of w, a, and o¢. 


13 On account of symmetry the shear stresses in the coordinate directions are of course zero. 
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We remark that the so-called membrane theory of axially symmetric shells results 
from the above theory when all non-linear terms in o», a4, u, and w are rejected. 

We have already stated in the introduction that the solutions of the variational 
equations (3.2) and (3.3) are unstable when the “radial” stress og is negative (i.e., 
when it is a compressive stress). On the other hand, it was stated that the solutions 
of the linear sheet theory are stable whether g¢ is positive or negative. The coriclusion 
regarding the instability in the non-linear case results immediately from the fact that 
the Legendre condition on the second variation of V is not satisfied if a» is negative, 
which means that the extremals do not render V a minimum in this case. The Le- 
gendre condition™ for a minimum in our case requires that the quantity A given by 


A = FupuP wzw; _ Faye (3. 7) 


should be positive. The quantity F is the integrand in (3.1) and subscripts denote 
differentiations. It turns out that the quantity A can be expressed in the form 
A=47°h?E(1—v*)—'x’o9. The right hand side has the sign of o¢, and consequently the 
Legendre condition is violated at all points where a» is negative. 

In the special case of the spherical sheet, it is possible to put the sign of g¢ in rela- 
tion to that of the applied pressure p. If the boundary conditions are specialized in 
such a way that the solution of the linear sheet theory results, we know [cf. (1.14) | 
that p and o¢ are opposite in sign, so that the solutions in this case are unstable when p 
is positive, i.e., when p is directed toward the center of the sphere. Since it is not pos- 
sible to give the solutions explicitly in the general non-linear case, we have not been 
able to prove readily that o» and p are opposite in sign in this case; but if the dis- 
placements remain small there can be little doubt that p and o¢ differ in sign in these 
cases also. In our further discussion of the spherical sheet we have therefore assumed 
always that p is negative, i.e., that it is directed outward from the center of the 
sphere. 

The linearized sheet theory results from (3.1) when all terms of degree higher than 
the second in u and w and their derivatives are neglected at the outset. If this is done, 
the Legendre condition for the resulting variational problem becomes A= Fy,u,>0, 
with 

2axEh 
A = ——— (3.8) 


1 — p? 


which is always positive, since x (the coordinate measuring the distance from the axis 
of the sheet) is always positive. Hence the Legendre condition is always satisfied in 
the case of the linear sheet theory, and we expect all solutions to be stable. The reason 
for the stable character of all solutions given by the linear theory, as contrasted with 
the unstable character of some of the solutions given by the non-linear theory, is that 
the linearization is equivalent to the imposition of a constraint powerful enough to 
cause stability in all cases. 

4. The spherical sheet. In the special case of the sphere we may write = Ré, 
p=R, and x=R sin 0, R being the radius of the sphere. The differential equations for 
the sphere are 





4 See, for example, R. Courant and D. Hilbert, Methoden der mathematischen Physik, vol. 1, Julius 
Springer, Berlin, 1931, p. 184. 
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= (os sin) = —-(— +1) on sin 0 + 04 c08@ (4.1) 
- og sin 0) = els u) oe sin os cos 8, ‘ 

1 d dw Rp , 

—- | on(= + w) sin | = - (+ + 0+ 04) sin 0. (4.2) 
R dé d0 h 


The system of equations is completed by the two strain-displacement relations (2.7) 
and (2.8) and the stress-strain relations (3.4). 

As boundary conditions at a fixed edge 6=@» we have [cf. the remarks preceding 
(3.5) and (3.6) | 


w= 0, (4.3) os — vog = 0. (4.4) 


In case the sheet has no hole at the axis, we require the solutions to be regular at 0=0. 
It is of interest to consider the special case of the complete sphere, in which the 
boundary conditions would become regularity conditions for 6=0 and @=7. In case 
we assume that the load is constant, it is readily verified that a solution of our dif- 
ferential equations which satisfies the regularity conditions is u =0, w= pR?(1—v)/2Eh, 
o9=04= —pR/2h. It could also be shown that this is the unique solution to this prob- 
lem. We observe that this solution is identical with that furnished by the linear sheet 
theory (a) of shells. In other words, the non-linear terms have no influence on the 
solutions for the full sphere in case the applied pressure p is constant. If p is not con- 
stant, however, the non-linear terms will influence the results for the full sphere. 

5. The spherical segment with small curvature. The differential equations of the 
Féppl-Hencky theory for the deflection of a radially symmetric plane sheet are con- 
tained in the above equations as a limit case. We need only allow R to tend to infinity 
while @ approaches zero in such a way that the product R sin 6 approaches a finite 
limit 7, and r is thus the polar coordinate which locates points in the plane sheet. 
The resulting equations (1.1)’ and (1.2)’ have already been given in the introduction. 
In passing to the limit, one observes particularly that the term u in the second paren- 
thesis on the right hand side of (2.7) drops out, so that the non-linear term reduces to 
4(dw/dr)*. As a consequence of this, the variational equations for the case of the plane 
sheet are much simpler than (4.1) and (4.2), since the terms corresponding to the 
first term on the right hand side of (4.1) and the term wu in the parenthesis on the left 
hand side of (4.2) disappear. 

It is clear that we could also simplify our equations for the spherical sheet quite 
considerably by omitting the non-linear terms involving u in the expression (2.7) 
for €9. It would seem fair to expect that such a simplification would be justified for the 
special case of rather flat spherical sheets. We recall that the choice of the expressions 
for the strains in terms of the displacements was in any case somewhat arbitrary. At 
the beginning, we might have considered the displacement u as negligible compared 
with the quantity dw/d@, since we expect the order of magnitude of the displacement 
w in the direction of the applied load to differ from that of the displacement u. In 
other words, it may well be that this term could be neglected even for sheets of rather 
large curvature.” In what follows we shall neglect this term. 





6 It might be noted that the limit problem which leads to the boundary layer phenomena (to be 
treated in the next section) is the same whether the terms in the displacement # under discussion here are 
retained or not. This is another valid reason for considering these terms to be negligible in most cases. 
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One result of this assumption is that the variational equations no longer contain 
the function u, but only w, og, and og. We can obtain a third equation in these same 
quantities—a “compatibility equation”—by eliminating u from (2.7) and (2.8) and 
then replacing €, and ¢, by a» and a, through use of the stress-strain relations. The 


result is the system of equations 


d 

















— (o,sin #0) = co, cos 8, (5.1) 
c 
1 =| dw . | |" + ‘ | me (5.2) 
— —| oo — sind | = — | — + 09 + og | sin 8, $.2 
R dol a0 nee 
d y dw 1 /dw\? 
—— (o, tan 0) + (1 + » tan? d)op = —]| w tan? 6 + tan 6 — + —|[—-] |]. (5.3) 
d@ R dé 2R \ dé 
0.5 
og] POWER SERIES 
k==1.56 x 1073 \uimit SOLUTION ,k=0 
XN 
- BOUNDARY LAYER \. 
w SOLUTION K==1.56x103\\, 
kR 
02 \ 
\ 
\ 
\ 
\ 
‘ \ 
\ 
\\ 
0.0 , 
0.00 004 0.08 0.12 0.16 0.20 


& (IN RADIANS) 
Fic. 3. Normal displacement. 


These equations are identical with Eqs. (1.11) which served as the basis for the dis- 
cussion of the curved sheet theory in the introduction. 

We are interested in solving the differential equations (5.1), (5.2), (5.3) for the 
case of a spherical segment without a hole about the axis @=0 and with a fixed edge 
at 8=6 . This means that we require the solution to be regular at 6=0 and to satisfy 
at the edge 6=6@, the conditions 


w= 0, (5.4) o4 — vog = 0. (5.5) 


6. Numerical solution of the boundary value problem for the flat spherical seg- 
ment. Our principal object in this paper is to present a new theory of thin sheets and 
to compare and contrast it with other theories, rather than to give numerical solu- 
tions for the resulting boundary value problems. However, we have obtained approxi- 


1945] NON-LINEAR THEORY OF CURVED ELASTIC SHEETS 261 


mate numerical solutions of the boundary value problem formulated in Eqs. (5.1) to 
(5.5) of the preceding section, and will report the results briefly in this section. 
The graphs of Figs. 3, 4, and 5 indicate the results of an approximate solution™ 


LiIMiT SOLUTION, k=0 
UTION 






























5 
POWER SERIES, > 
a ~ 
-k=1.56x1073 ee hy 
0.44 
BOUNDARY LAYER 
SOLUTION, -Kk =1.56 x 1073 
0.3 
-& 
kE 
0.24 
\ 
\ 
\ 
0.14 \ 
0.0 3 = . J . 
0.00 0.04 0.08 0.12 0.16 0.20 
J (IN RADIANS) 
Fic. 4. Circumferential stress. 
0.80 BOUNDARY LAYER SOLUTION,-k =1.56 x 1073 
Seeabetiaenenenseres } 
re. LIMIT SOLUTION, k=0 
~~ 
POWER SERIES——~ ~~ 
-k =1.56 x1073 ~ 
0.4874 \ 
\ 
\ 
\ 
\ 
\ 
0.46; \ 
Og \ 
\ 
KE \ 
0.444 \ 
\ 
\ 
\ 
0.424 
0.40 . _—e . . 
0.00 0.04 0.08 0.12 0.16 0.20 


9 (IN RADIANS) 


Fic. 5. Radial stress. 


16 The solution was obtained in the form of a development in powers of the independent variable @. 
Only four terms in the series (which appear not to converge very rapidly) were retained in calculating 
coefficients. On the graphs these solutions are marked “power series, — x = 1.56 X 10~*.” 
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for the case in which the value 6) of @ at the edge of the sheet is 0.2 rad., and the di- 
mensionless parameter!’ x= ~R/Eh has the value —1.56X10-*, and v=0.3. (This was 
one of the cases treated as part of the project mentioned in the footnote at the be- 
ginning of the paper.) The strain everywhere in the sheet is about —}(1—v)« in 
value. The graphs show values of w/xR, —o,/xE, and —oo/xE as functions of 6. 
Each graph contains three curves: a dotted curve giving the result of our approxi- 
mate solution; and two solid curves which refer to the limit solution obtained as x—0. 
The method of obtaining the limit solutions is explained in the next section. One ob- 
serves that the curves marked “boundary layer solution” approximate those of our 
numerical solution rather well, at least for the stresses. We have some reason to think 
that a more accurate solution of the boundary value problem would show the bound- 
ary layer solution to be a better approximation to the actual solution than our graphs 
indicate. We note that the curves marked “limit solution, x=0” are those which 
would be obtained from the linear sheet theory. 

7. The boundary layer problem. A boundary layer effect has already been men- 
tioned a number of times in connection with our boundary value problem. The graphs 
of the solutions in the preceding section furnish a hint regarding the character of this 
phenomenon. The solutions in the interior portion of the sheet appear to be relatively 
constant, approaching there the values furnished by the linearized theory (i.e., those 
of the theory usually called the membrane theory of shells). However, toward the 
edge of the sheet, the solutions appear to change rather abruptly. This is consistent 
with the repeatedly mentioned fact that the condition for a fixed edge cannot be satis- 
fied in the linearized theory. The purpose of the present section is to treat this bound- 
ary layer effect explicitly. 

A necessary step in any treatment of boundary layer phenomena’ consists in the 
introduction of appropriate new variables and parameters. In the present case it is 
convenient to introduce new dimensionless dependent variables replacing €, €4, 76, 7, 
w and u by the relations 


59 = oo/ Ex, S¢ = o4/ Ex, (7.1) 
€9 = €9/k, €y = €4/k, (7.2) 
w = w/Rx, wa = u/Rx, (7.3) 


in which the important dimensionless parameter x is defined by the relation 
xk = pR/Eh. (7.4) 
We assume here that the applied pressure p is constant. In terms of the new quanti- 


ties, the fundamental differential equations (4.1), (4.2), (2.7), and (2.8) become, in 
order 


= (sg sin 6) = (+n) se sino +5, cos 8, (7.5) 
dé dé 
d dw ; 

| o(S +4) sino | = - [1 + so + sy] sin 4, (7.6) 


17 In the next section it will be seen that 0) and « are the only essential parameters, once the value of 


the Poisson ratio »v is fixed. 
18 Such boundary layer effects have been well-known for many years in fluid mechanics. They occur 


also in problems in elasticity other than those considered in this paper. (See, for example [7, 8].) 
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(E-es(es) : 
a) tee dae 


es = (u cot @ — w). (7.8) 


€e 


To obtain a complete system of equations we add the stress-strain relations: 
€o = So — Sg, (7.9) Cy = Sy — Sp. (7.10) 


As boundary conditions we require all quantities to be regular at @=0, while at 6=4) 
the condition of a fixed edge is prescribed, 


w = 0, (7.11) so — vy = 0. (7.12) 


We now observe that if x is allowed to approach zero in these differential equa- 
tions, the result is a set of differential equations for the limit quantities which are 
identically the same as those of the linear sheet theory’® (when formulated in terms 
of our dimensionless variables), 

d 
— (spsin@) = s,cos8, (7.13) Sots, = — 1, (7.14) 


dé 
ee = — — a, (7.15) és = uw cotd — w. (7.16) 


Obviously, the boundary conditions (7.11) and (7.12) cannot be imposed in this limit 
problem. In fact, the solutions of (7.13) and (7.14) are completely determined by the 
regularity conditions at @=0 alone. This solution is, as we know, se=s,=—}, 
w=3(1—v), »=0. In the limit, therefore, the boundary conditions at the edge, in 
general, will not be satisfied. It follows that the solution of the boundary value prob- 
lem formulated in (7.5) to (7.12) will not converge uniformly at the boundary to the 
solution of the limit problem as x-0, and this is the essential characteristic of a 
boundary layer effect. 

It is, however, reasonable to expect that the solutions do converge in the interior 
(i.e., for OSO0S0¢ <@o, where 6¢ is a constant) as x-0 to the solutions of the limit 
problem for x=0. The graphs of the preceding section confirm this to some extent. 

It is possible to give an explicit treatment of the boundary layer effect. Such a treat- 
ment can be obtained through the introduction of a new independent variable which 
replaces @ and which depends on «x in such a way that the solutions are made to con- 
verge uniformly at the boundary in the limit as x0. What one wants, roughly speak- 
ing, is to stretch the boundary layer as x0 in such a way that its width does not 
shrink to zero. In our case, this can be accomplished by introducing as a new inde- 
pendent variable the quantity 7 defined by the relation 


(@ — 6). (7.17) 





i 
V-k 





19 Tt is perhaps of interest to observe that the limit situation characterized by x0 can be achieved 
by allowing the pressure p to approach zero. However, if p is simply set equal to zero in the original dif- 
ferential equations (4.1) and (4.2), the order of the system is not decreased. Thus the introduction of 
new dependent variables through division of the original ones by « is an essential step in the treatment 
of the boundary layer effect. 
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The reason for the minus sign under the radical is that we wish to consider only cases 
for which the solutions are stable, which means cases in which x is negative. (Cf. the 
remarks at the end of Section 3.) 

If we introduce the new independent variable in Eqs. (7.5) to (7.8) and then allow 
x to tend to zero, we obtain the set of limit differential equations 


dso A . d dw ; 

ak (7.18) bs WS) = (1+se+s,), (7.19) 

dn dn dy 

; du 

0= : ’ (7.20) €y = UCOt Op — w = Sg — ve, (7.21) 
dy 


for the range —:0© <m<0, where y7=0 corresponds to the edge of the sheet. This sys- 
tem of equations, which has the same order as the original system, yields the bound- 
ary layer “resolution” which we seek. The boundary conditions at 7=0 are given by 
(7.11) and (7.12). At y= — © we expect all quantities to tend to the values furnished 
by the solution of the interior limit problem given above. Thus we expect w to ap- 
proach the value $(1—v) as n> — &. 

Since the boundary layer differential equations have constant coefficients, they 
are readily solved by exponentials. One finds, for example, that w satisfies the differ- 


ential equation 





d*w 
-— 2 = —1+»,, (7.22) 
dn? 
so that the homogeneous equation is solved by real exponentials.?® The solution of 
(7.22) which satisfies the conditions at 7=0 and n= — & is 
w = 4(1 — v)(ev2" — 1). (7.23) 
The results for the other quantities are easily found to be 
zp =0, (7.24) So = — j, (7.25) 
so = —4{1+ (1 — v)(e?"— 1)}. (7.26) 
The graphs of Figs. 3, 4, and 5 contain in each case a curve marked “boundary 
layer solution, x= —1.56X10-*.” These curves were obtained from (7.23), (7.25), and 


(7.26) by reintroduction of @ as a variable through use of (7.17) with x= —1.56 X10-%. 
Comparison with the curves for the numerical solution of the original boundary value 
problem indicates that such a “compressed” boundary layer solution may furnish a 
fairly good approximation to the values of w and o» near the edge of the sheet if x is 
not too large. 
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—NOTES— 


A FORMULA FOR THE SOLUTION OF AN 
ARBITRARY ANALYTIC EQUATION* 


By D. R. BLASKETT anp H. SCHWERDTFEGER (University of Adelaide) 


In this note a proof will be given of a formula which has been stated without 
proof by E. Schréder.' The fact that he has expressed some doubts as to its general 
validity may have caused it to fall into oblivion although it seems to be of some theo- 
retical and practical interest, like several other results of Schréder’s to be mentioned 
below. Since moreover the true nature and simplicity of the formula is rather con- 
cealed in Schréder’s discussion, it may be worth while to enter into the matter again. 

The formula in question is a consequence of the following theorem which at once 
describes its exact realm of validity: 


THEOREM I. Let w=f(z) be an analytic function, regular in a domain A of the com- 
plex z-plane. Let a be an interior point of A and a simple root of f(z); 


f(a) = 0, f(a) # 0. 


Further let zo be a point in A “not too far from” a (in practice 201s a first rough approxt- 





mation to the root a). Then, denoting by 2=f~'(w) the inverse function of f(z) one has 
= f(20)’ (@’f-"(w) df~'(w) 
eB LX (LO) ee (- se) 
p= y! dw” w=f(z9) .« dw w=f (29) 


where the exponential function operates symbolically on the differential symbol. 


In this form the theorem is a corollary to the main theorem on the analyticity 
of the inverse of an analytic function. To prove (1) we make use of the fact that be- 
cause f’(a)0 the function f(z) is simple (schlicht) in a certain neighborhood of 
the point z=a.? Hence the inverse f~'(w) exists in a circle K of radius p>0 round the 
point w=0. Moreover it is analytic in K. Thus we may choose 2» in a neighborhood 
of a such that a circle k round wo=f(zo) (in the w-plane) contains 0 and is wholly 
contained in K. In k the function f-'(w) is given by its power series 


1 (CL) “(w — f(zo))” 
w=f (Z0) 


Hd = 2. — 
yao V! 
whence for w=0 follows the formula (1). Evidently, for zo, any point z near @ can 
be chosen for which w=f(z) lies in a circle of radius p/2 round 0. 
To obtain Schréder’s formula we introduce the operators 6“ (u=0, 1, 2,-- - ) 
defined as follows: 





dw” 


* Received Oct. 30, 1944. 

1 E, Schrider, Uber unendlich viele Algorithmen zur Auflésung der Gleichungen, Mathematische An- 
nalen 2, 317-365 (1870); in particular cf. §§4 and 5 of this paper. 

2 Cf. G. Julia, Principes géométriques d’analyse, Premiére Partie, Paris 1930, p. 16-17. 
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5°f(z) = 5'f(z) . = ( ). 5*f(z) 
ne f(z) = "f(z) = 
f'(2) f'(2) f'(2) 


Siffe\ =— d*f-\(w) J = i a 
ose) -( dw? his *fe) -( dw’*} esi 


which is easily shown by induction. Therefore the expansion (1) is identical with 


f'@) 


Then 











a = 20 + u (- iy — Ks wl —— (P- f(z om zo 
1 fe) f’'(z0) , f(zo)® f’(20)f’”’ (zo) — 3f’"(z0)? 
= i6 50 2 + + a rae 2 
- f’ (20) 2! f'(z0)* 3! f'(%0)* ” 


This was not the way in which the formula (2) was actually discovered. It ap- 
peared rather as a plausible consequence when dealing with another method of ap- 
proximate solution of an analytic equation. Although no proof could be given on 
the basis of this method, it turns out to be of some interest here as it shows that the 
partial sums of the infinite series (1) or (2), i.e., the expressions 





®, (2) -X(- 1) LO” ar) - > (—1 yp Le” A yl 3) 
! dw’ a 
can be used for an iterative approximation of the root a; 
THEOREM II. For each n=1,2,-- + the recurring sequence 
ao = Zo, Om = P,(am—1), (m =1,2,---) (4) 


has the root a as its limit and then also 
&,(a) = a, &/(a) = 0,---, & (a) = 0. (4’) 


This is one of the results of Schréder who has discussed in detail (I.c.) several 
such iterative algorithms (“Algorithmen erster Art”) and has thus obtained formulae 
of considerable practical interest. We propose here another treatment of the problem 
of iterative approximation which leads immediately to a proof of Theorem II. 

We make use of the fact that if a is a simple root of the equation f(z) =0, a func- 
tion ¢(z) can be found for which a is an attractive fixed point, i.e.,‘ 


g(a)=a, |¢'(a)| <1. (5) 


Then the recurring sequence ao, a1 =¢(ao), 22 =o(a1), - - - , if ao is not too far from a, 
is convergent and has @ as its limit. To strengthen the convergence we may replace 
the inequality in (5) by the condition ¢’(@) =0. Such a function ¢(z), involving an ar- 
bitrary funct tion h(z), is for instance 








3 In pe articular the reader’s attention may be called to the formulae (A> ) and (BY) i in §12 (p. 352) of 
Schréder’s paper. These formulae are very expedient for the computation of mth roots if a high degree of 
accuracy (e.g. more than 20 correct figures) is required. A similar algorithm has recently been given by 
V. A. Bailey in a brief expository article Prodigious calculation, Australian Journal of Science, 3, 78-80 
(1941), by which our attention has been drawn to the present subject. 

4 Cf. G. Julia, loc. cit., p. 23. 
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f(a) h(z) — i 
2) = 2 + f(z) ——— 
(2) f(z) ye) 


Thus we may impose more rigid conditions further strengthening the convergence, 
viz., 


¢'(a) = 0, g(a) = 0,--+ , o™(a) = 0, (6) 


where 7 is any positive integer.® 
A function ¢(z) satisfying the conditions (5) and (6) can be obtained in the follow- 
The conditions (6) will be satisfied if the derivative of ¢(z) appears in the 


¢'(z) = (f(z) "g(2)f'(2), 
the undetermined function g(z) being regular at z=a. It remains to adapt g(z) to the 


condition ¢(a@) =a. One has 
(2) = f e@yreoy' eas = f wg F(w)aw 


whence, by repeated integration by parts, it follows that 


1)” 
(z) = n! a. —— (f(z))*~"go41(s) 


ym eT 


ing way. 
form 


where g,(z) is the u-fold iterated indefinite integral of g(f-!(w)) for w=f(z). Thus, for 
z=a one has 
o(a) = (— 1) "1! gn4i1(a). 
Therefore 
(= a)" 


n! 


En + 1(z) _ 


will give a function ¢(z) which has all the desired properties. In this way one obtains 
the function ®,(z) of (3), and it is evident that this function has the properties stated 


in Theorem II. 





5 From a letter of Professor V. A. Bailey we have learnt (in May 1941) that this problem has been 
dealt with in some special cases by E. Netto in his Vorlesungen tiber Algebra vol. 1, Teubner, Leipzig, 1896, 
p. 300. In the same letter Bailey has given an elegant solution of the problem which, however, does not 
suit our present purpose. Further he has drawn our attention to the paper by L. Sancery, De la méthode 
des substitutions successives pour le calcul des racines des équations, Nouvelles Annales d. Math. (2) 1, 
305-315 (1862), which, however, was not accessible to us. 


THE CAPACITY OF TWIN CABLE* 
By J. W. CRAGGS anp C. J. TRANTER (Military College of Science, Stoke-on-Trent, England) 


1. Introduction. The problem of determining the capacity of two long parallel 
cylindrical conductors can be easily solved by the use of a conformal transformation." 
A simple extension of the method gives the result for the case in which each conductor 


* Received April 16, 1945. 
1F, B. Pidduck, A treatise on electricity, Cambridge University Press, 1916, p. 77. 
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is surrounded by a dielectric sheath whose boundary is a member of the coaxial sys- 
tem of circles defined by the boundaries of the conductors. The case in which the 
sheaths are concentric with the conductors is of much greater practical importance 
and in many types of cable the sheaths are actually touching. In this paper we give 
the derivation of the potential distribution for this latter case together with a prac- 
tical method for the calculation of the capacity. 

2. Statement of the problem. We consider the symmetrical problem of two circu- 
lar wires each of radius R; surrounded by concentric touching sheaths of radius R: 
and dielectric constant K,;, the whole being immersed in an infinite medium of dielec- 
tric constant Ko. 

For infinitely long straight wires, the problem reduces to the determination of 
potentials Vi, V2 satisfying: (7) the equations 


V*Vi 


0, (1) 


for Ri Sr= Ro, and 
V°V2 = 0, (2) 


in the region between the circle r= R: and the line x =0, where V? is the two dimen- 
sional form of Laplace’s operator, the polar coordinates 7, 6 are based on the centre 
of one of the conductors and the cartesian coordinates x, y have origin at the point 


y 














of contact of the sheaths and axes as shown in Fig 1; (77) the boundary conditions 


Vi=1 (3) 

when r = R;, and 
Vi = Vo, (4) K,0V,/or = K,0V2/dr, (5) 

when r=Ro, and 
V2 = 0 (6) 


when x =0. Condition (6) is a result of the symmetry of the problem provided that the 
potential on the left hand conductor is taken as —1. 
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The capacity C per unit length of wire is then given by C=3Q, where Q is the 
charge per unit length on either conductor. 

The capacity per unit length is unaltered if we replace R: by unity and R; by 
R;/R2(=a) and we shall do this in the subsequent work. 

3. The analytical solution. It is natural to express the potential V; in the polar 
coordinates defined above. We therefore write 


r ta) r n a n 
Vi=1+ Blog— + > {(<)- (<) os, cos n6, (7) 
a n=1 a r 


this being the most general solution of (1), symmetrical in 6 and satisfying (3). 
Conformal representation by the use of 





[—s (8) 
—_— 1n a -« E = 

z+iy 1+  re** 

transforms the region r>1, x>0 into 0<£<}3, the boundaries x =0, r=1 becoming 
£=0, =} respectively. The general solution of (2) satisfying (6) and the conditions 
of symmetry is 


V2 -f f(t) sinh 2& cos 2nédt. (9) 
0 


The constants B, b, of (7) and the function f(t) of (9) are now to be determined 
from the boundary conditions (4) and (5). 
Now on r=1 (£=3) the relation (8) gives 








n = 3 tan 30 = 36 . (10) 
say, and 
OV : ‘~ OV 11 + BY) OV (11) 
— = — i sec? 49 — = — 2) —_. 
ar ’ wa . dg 
Thus (4) and (5) become 
=) 1 _— aq?" oo 
1— Bloga+ b b, cos n0 = f f(t) sinh ¢ cos Biédt, (12) 
n=1 a” 0 
7 _ +e - 
KB+K be — nb, cos nO = — (1 + ef t{(t) cosh ¢ cos Bidt, (13) 
n=1 a” 0 


where K = K;,/K2. 
Multiplying (12) by cos n@ (n=0, 1, 2,---) and integrating with respect to 0 


from 0 to 7, we have 


1 ® re) oo 
1— Bloga= —f f f(t) sinh ¢ cos Bididé = f e~*f(t) sinh édt (14) 
T 0 0 0 


* = cos Bi 
f cos Bid@ = 2f dB = re 
0 o 1+ 8? 


since 





and 
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2 2 ® ) ) 
—_——}, = - f f f(t) sinh ¢ cos Bt cos nOdtd6 = f e~*f(t) sinh ¢7,(@dt, (15) 
TT 0 0 0 




















where 
2e! . 4e' © cos Bt 
T,(t) = -f cos n6 cos Bid8 = —— cos (2n tan! B)dB 
T 0 Tv 0 + B? 
> —, 
= ( — 1)? ———— (24)""”. (16 
p=0 (n — p)! 
Applying Fourier’s integral theorem to (13), we obtain 
t 4B ° cos Bt 2 1+ a™* cos Bt 
— — f(t) cosh ¢ = — B ee nb,, i cos n@ dp 
K x Jo i+ ps" T nal 1 + 6? 
Cs) 1 4+ a?" 
= 2Be-* + e* >, ——— nb, I,(2). (17) 
n=1 a” 
Equations (15) and (17) lead to 
1 — a? 2 1+a™* 
— —— b, = Ba, + ), ——— nbAuy, (18) 
Ka? n=1 a” 
where 
> tanh ¢ 
a, = 2f e~?t I ,(é)dt, (19) 
0 
and 
4 tanh ¢ 
Anp = f ” is ; T,(t)T,()dt. (20) 
0 
Finally (14) and (17) give 
» tes is i+e 
(B log a — 1) = 2B f " os ; —_ a ae NbpQn. (21) 
0 


The infinite set of equations (18) gives the values of the coefficients b in terms of B. 
Substitution in (21) yields an equation for B and the capacity per unit length C can 


then be determined from C= —1K3,B, since 
Ky al OV; 
2c =Q=-— —— adé. 
4n 0 or 


This completes the analytical solution. 

4. Method of computation. In practice, a good approximation to the capacity 

may be obtained by retaining only a finite number m of the coefficients b,. Eliminating 

the m quantities [(1+a’")/a"]nb, (n=1, 2, - - - , m), from the (m+1) equations (18) 
ar.d (21), and writing 

1 — a™* 


= ————————, (22) 
nK(1 + a?) 
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we obtain 





| Aa +4 ; ere Bus ay 
Ag Asoo + Y¥2°°° Aom Qe 
=0, (23) 
A m1 Ame*** Amm + Ym Qm 
a vy 
a1 Qo: Om > (= — log a) + 4 log = 


since 
2s tanh ¢ T 
f e-## ___— dt = log — - 
0 t z 


The quantities a,, An, can be easily computed from Eqs. (19), (20) and (16) with 
the help of the result? 


*” 1 
f (2t)"e-** tanh tdi = ni {(1 — a iG +i) - i} , 


where {(m) is the Riemann Zeta-function, tabulated for integral m in J. Edwards, 
“The integral calculus,” vol. 2, Macmillan, London, 1922, p. 144.* 

The solution for K =1 differs from the well-known exact solution for this case by 
less than 0.2 per cent, when only the first three of b, are retained, provided that aS}. 
For larger values of K and a it may be necessary to retain more terms to achieve the 
desired accuracy but for practical values the amount of computation required is not 
excessive. 

2 When n=0, the result reduces to fie? tanh tdt=log 2—}. 

3 A four-figure table is given in E. Jahnke and F. Emde, Tables of functions, Dover Publications, 


New York, 1943, p. 273. 


LARGE DEFLECTION OF CANTILEVER BEAMS* 
By K. E. BISSHOPP anp D. C. DRUCKER (Armour Research Foundation) 


The solution for large deflection of a cantilever beam! cannot be obtained from 
elementary beam theory since the basic assumptions are no longer valid. Specifically, 
the elementary theory neglects the square of the first derivative in the curvature 
formula and provides no correction for the shortening of the moment arm as the 
loaded end of the beam deflects. For large finite loads, it gives deflections greater than 
the length of the beam! The square of the first derivative and correction factors for 
the shortening of the moment arm become the major contribution to the solution of 





* Received April 6, 1945. 
1 This problem was considered by H. J. Barten, “On the Deflection of a Cantilever Beam, 


of Applied Math., 2, 168-171 (1944). Previously an approximate solution had been obtained by Gross und 
Lehr in Die Federn, Berlin VDI Verlag, 1938. 


” 


Quarterly 
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large deflection problems. The following theory which utilizes these corrections is in 
agreement with experimental observations. 

The derivation is based on the fundamental Bernoulli-Euler theorem which states 
that the curvature is proportional to the bending moment. It is assumed also that 
bending does not alter the length of the beam. 

Considering a long, thin cantilever leaf spring, let L denote the length of beam, 
A the horizontal component of the displacement of the loaded end of the beam, 6 the 
corresponding vertical displacement, P the concentrated vertical load at the free end, 
B the flexural rigidity, that is B=EJ, when cross-sectional dimensions are of the 








4 
rs 
Nu 


Fic. 1. 


same order of magnitude, and B=EI/(1—v?) for “wide” beams, where » is the Pois- 
son ratio. The exact expression for the curvature of the elastic line may be stated 
conveniently in terms of arc length and slope angle denoted by s and ¢, respectively, 
so that if x is the horizontal coordinate measured from the fixed end of the beam, the 


product of B and the curvature of the beam equals the bending moment M: 
do 

B— = 
ds 


P(iL—x-A)=M (1) 


or 
d*$ P dx P 4 (2) 
—_ = — — — = — — 0s @, 
ds? B ds B 
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whence 


oe Roe oy 
iw 68 ~=s 


(3) 


The constant C can be evaluated directly by observing that the curvature at the 
loaded end is zero. Then if ¢o is the corresponding angle of slope 


dp calf in g) 12 (4) 
7s = 3B sin do sin @)*!*. 


The value of ¢o cannot be found directly from this equation but it is implied by the 
requirement that the beam be inextensible, so that 


oP pt $0 _/ PL2\1/2 
— ds = f ‘sin do — Si “1/2dg = 2( ) F 5 
/=f (sin @) — sin @)~'/"de = y - (5) 


In order to evaluate this elliptic integral, denote PL?/B by a? and let 





1 + sin @ = 2k? sin? 6 = (1 + sin go) sin? 0 
Then 


x/2 
a -f (1 — k?* sin? 6)—'/d6, 
61 


The next step is to represent the deflection 6 in terms of @ and an elliptic integral 
Since 


(6) 


sin 0, = —/2/2k. (7) 


dydp dy ; 
—-—— + —— = §i @, 
do ds_ ds 
and since we have d¢/ds from Eq. (4), 
dy /2P ng)? = si 
—A/ — (sin do — sin 2 = sin ¢. 
ie B do $ ¢ 
Thus 


y "B fe sin od 
6 -f dy -4/—f wt 
0 2F 


0 (sin do — sin o)*/? ; 





With the aid of Eq. (6) we obtain 


6 V 3 | $0 sin odd * (2k? sin® 6 — sin? 6 — 1)d0 
Lady (sings — sing)? 





61 a _ — & sin? sin? 6)}/? 
This equation can be split up into complete and incomplete elliptic integrals of the 
first and second kinds. In the notation of Jahnke and Emde, 


6 


 to- [F(k) — F(k, 0:) — 2E(k) + 2E(k, 61) | (8) 
a = F(k) — F(R, 4), 
so that 
é 2 
— = 1—— [E(k) — E(k, 6;)]. (9) 
a 


1945] H. J. BARTEN 275 


The horizontal displacement of the loaded end is calculated from Eqs. (1) and (4) with 
x=0 when ¢=0. Thus 


P(L — A) 5(2) By/™ oi ye 
. ds ied B ute 


L-A v2... 
7 = one: (10) 


or 





From Eq. (6) we have sin ¢9=2k?—1. 

Numerical results can be obtained by: (1) selecting values of k corresponding to 
tabulated values of the modular angle in the elliptic function tables and (2) determin- 
ing 6; and a from Eq. (7). After this has been done, 6/Z and (L —A)/L can be calculated 
from Eqs. (9) and (10) and plotted against a? =PL?/B. The results of these calcula- 


tions are shown in Fig. 1. 


CORRECTIONS TO MY PAPER 
ON THE DEFLECTION OF A CANTILEVER BEAM* 
QUARTERLY OF APPLIED MATHEMATICS, 2, 168-171 (1944) 


By H. J. BARTEN 


This paper is correct up to the equation 


L 
Or -f as cos 6 ds. 
0 


d6 ;, 
— = aL cos 6, 
dL 


is incorrect since @ is not only a function of L, but is also a function of s. This error 


makes Eqs. (9), (11), and (12) incorrect. 
Using the relation 


The next step 


dé 
— 2 Gist — © 

ds 
and the various steps used in the original paper, we find that 


aL = F( 2, =) — F(k, 8). 


By using 6 as an independent variable we can calculate corresponding values of k and 











* Received June 25, 1945. 








aL? which in turn are used to 
find corresponding values of 
F, and F,. The corrected 
curves thus derived are shown 
in Fig. 1. 

The author wishes to 
thank M. M. Johnson of 
Washington, D. C. and D. 
C. Drucker of the Armour 
Research Foundation for 
pointing out these discrepan- 
cies. 
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